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Abstract 

The  United  States  (US)  military  and  its  allies  have  found  themselves  hghting 
in  asymmetrical  wars  against  terrorist  and  insurgent  groups.  Allied  forces  have  been 
hesitant  to  attack  targets  at  certain  locations  due  to  the  probability  of  injury  or  death 
to  innocent  civilians.  Low  fragmentation  weapons  with  concentrated  pressures  that 
dissipate  quickly  and  the  addition  of  aluminum  {At)  particulates  in  explosive  mixtures 
will  provide  the  capabilities  needed  to  develop  weapons  needed  to  fight  current  and 
future  asymmetric  warfare. 

The  issue  with  such  weapons  is  that  the  physics  of  multi-phase  explosives  are 
not  well  understood.  The  inclusion  of  metal  particulates  in  explosives  represents  a 
particularly  computationally  challenging  physical  environment,  in  that  such  flows  are 
“two-phase”  flows.  In  “two-phase”  flows,  gas  physics  and  solid  particle  physics  are 
simultaneously  simulated.  Since  the  two  phases  represent  different  states  of  matter, 
matching  of  physical  conditions  for  the  two  different  phases  is  difficult  to  accomplish 
computationally.  The  assumptions  made  for  the  development  of  a  computational 
code  that  simulates  multi-phase  explosives  are  presented  in  this  document. 

This  research  focused  on  analyzing  the  effects  of  aluminum  in  high  explosives 
such  as  PBX9501  using  a  computational  simulation  code.  The  Multi-Phase  Explosive 
Simulation  (MPEXS)  hydrocode  is  used  to  perform  simulations  on  a  wide  range  of 
metalized  explosive  cases,  where  the  varying  parameters  are  the  size  of  the  aluminum 
particles  from  9  to  34  microns  and  the  mass  fraction  of  aluminum  from  10  to  20 
percent.  The  mass  fraction  directly  correlates  to  the  amount  of  aluminum  present  in 
the  mixture.  The  MPEXS  code  was  evaluated  by  performing  several  studies  including 
a  data  convergence  test  to  determine  the  stability  and  limitations  of  the  code. 
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The  important  explosive  detonation  parameters  that  were  considered  to  evalu¬ 
ate  the  effects  of  aluminum  were  pressure,  density,  detonation  velocity,  and  Run-to- 
Detonation  (Run2Det)  time  and  distance.  The  results  showed  that  the  gas  pressure 
steadily  increased  as  the  A1  particle  diameter  increased.  When  the  A1  mass  fraction 
was  increased,  the  gas  pressure  for  each  of  the  particles  decreased.  The  gas  density 
displayed  the  same  trends  as  pressure.  Density  steadily  increased  as  the  A1  particle 
size  increased,  and  density  decreased  when  A1  mass  fraction  increased.  Detonation 
velocity  steadily  decreased  as  the  A1  particle  diameter  was  increased,  and  detonation 
velocity  decreased  as  the  A1  mass  fraction  increased.  It  was  concluded  that  varying 
the  aluminum  mass  fraction  had  a  larger  effect  on  Run2Det  properties  than  varying 
the  aluminum  particle  diameter  for  both  Run2Det  distance  and  time.  In  general, 
increasing  the  A1  mass  fraction  will  decrease  the  the  time  and  distance  for  the  deto¬ 
nation  to  reach  steady-state. 

The  addition  of  metal  particulates  in  explosive  mixtures  increases  the  density  of 
the  shock  wave,  causing  a  higher  pressure  in  the  shock.  The  high  pressure  in  the 
shock  is  devastating  and  will  incapacitate  adversaries  in  physical  proximity  of  the 
detonation  of  the  munition.  The  concentrated  pressures  will  dissipate  quickly  due  to 
the  high  density,  providing  a  short  proximity  detonation.  This  research  signihcantly 
contributes  to  both  current  explosive  simulation  analysis  and  development  of  future 
explosive  formulations  for  better  munitions  that  will  save  many  lives. 
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SIMULATION  OF  METAL  PARTICULATES  IN  HIGH  ENERGETIC 

MATERIALS 

I.  Introduction 

Asymmetrical  warfare  has  presented  many  issues  to  the  US  and  its  allies  in  the  way 
wars  are  being  fought.  Allied  forces  have  been  hesitant  to  attack  targets  at  certain 
locations  due  to  the  probability  of  injury  or  death  to  innocent  civilians.  One  solution 
to  this  problem  is  the  development  of  more  advanced  weapons.  Low  fragmentation 
weapons  with  concentrated  pressures  that  dissipate  quickly  might  provide  a  solution 
that  eliminates  targets  with  less  or  no  collateral  damage.  The  addition  of  Al  particu¬ 
lates  in  explosive  mixtures  might  provide  the  capabilities  to  develop  weapons  needed 
to  hght  current  and  future  asymmetric  warfare.  Explosive  simulation  hydrocodes 
such  as  the  MPEXS  code  will  allow  researchers  to  study  multi-phase  explosives  and 
formulate  new  explosive  material  for  use  in  future  advanced  munitions.  Section  1.1 
provides  a  background  on  the  current  issues  of  asymmetrical  wars  the  US  is  facing. 
Section  1.2  introduces  the  topic  of  metalized  explosives,  and  discusses  the  assump¬ 
tions  and  limitations  of  computational  modeling  of  explosives.  Section  1.3  presents 
the  research  questions  that  were  explored  in  this  research.  Finally,  Section  1.4  gives 
an  overview  of  subsequent  chapters  in  this  document. 
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1.1  Background 


1.1.1  Asymmetrical  Warfare 

In  the  recent  Iraq  and  Afghanistan  wars,  the  US  military  and  its  allies  have 
found  themselves  hghting  in  asymmetrical  wars  against  terrorist  and  insurgent  groups. 
These  groups  have  resorted  to  strategies  and  tactics  that  involve  hiding  in  and  attack¬ 
ing  from  heavily  populated  zones,  hospitals,  schools,  mosques,  and  other  locations. 
This  type  of  asymmetrical  warfare  makes  allied  forces  hesitant  to  attack  such  locations 
using  conventional  weapons  due  to  the  risk  of  damage  to  non-combatant  targets  or 
threat  to  innocent  civilians.  Conventional  weapons  such  as  air-to-surface  missiles  and 
bombs  are  designed  to  destroy  targets  effectively  by  destroying  entire  buildings  and 
town-blocks  where  adversaries  are  located.  The  problem  with  asymmetrical  warfare  is 
that  the  US  military  can  no  longer  destroy  targets  using  conventional  weapons.  The 
US  needs  the  capability  to  eliminate  targets  with  diminished  risk  of  killing  innocent 
civilians. 

War  strategy  is  dependent  on  many  factors  and  it  is  always  changing  to  meet 
the  needs  of  war.  One  thing  that  has  not  changed  much  is  that  wars  have  been 
fought  in  a  conventional  way  for  a  long  time.  Conventional  warfare  can  be  dehned 
as  a  confrontation  between  two  or  more  groups  with  well-dehned  confrontation  lines. 
These  lines  or  boundaries  can  be  a  location  on  a  held  of  battle,  or  the  border  line  of  a 
country,  or  a  boundary  surrounding  the  location  where  the  opposing  side  is  located. 
Either  both  sides  come  to  accordance  on  where  the  battle  will  be  fought  or  one  side 
attacks  the  other  side’s  territory.  In  both  cases  there  is  a  clear  distinction  between 
the  groups/enemies.  The  purpose  of  conventional  warfare  is  to  weaken  or  destroy  the 
enemy’s  military.  When  destroying  or  weakening  the  enemy  is  not  possible,  either 
side  might  resort  to  asymmetrical  warfare  to  exploit  the  other’s  weakness. 

The  US  and  its  allies  have  shown  superiority  over  enemies  and  terrorist  groups 
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when  it  comes  to  conventional  warfare.  Operation  Desert  Fox  was  an  example  of  an 
effective  attack  campaign  that  destroyed  the  enemy’s  military  infrastructure  in  a  short 
period  of  time  with  no  allied  losses  [10].  Operation  Desert  Fox  was  the  “code-name” 
for  the  bombing  of  Iraq  on  December  16  -  December  20,  1998  [10].  The  invasion  of 
Iraq  in  2003  was  another  example  of  a  successful  campaign  in  a  conventional  type 
of  war.  From  20  March  2003  -  1  May  2003  the  coalition  forces  invaded  Iraq  and 
overthrew  the  Ba’athist  Iraqi  government  led  by  Saddam  Hussein.  The  campaign 
ended  when  Coalition  forces  captured  the  city  of  Baghdad.  To  put  in  perspective  the 
superiority  of  the  US  in  conventional  warfare,  coalition  forces  attacked  an  enemy  of 
equal  troop  size  and  defeated  the  Ba’athist  Iraqi  government  while  only  suffering  a 
few  hundred  losses. 

Terrorist  groups  and  Taliban  insurgents  learned  that  their  chances  of  defeating 
the  US  in  conventional  warfare  are  close  to  none  and  therefore  they  have  resorted 
to  hghting  with  asymmetrical  tactics.  These  adversary  groups  have  found  that  the 
US  will  not  attack  hospitals,  schools,  mosques  and  populated  zones;  therefore  these 
groups  use  such  locations  as  hiding  shelters  and  weapon  storage.  Asymmetrical  war¬ 
fare  presents  many  problems  for  allied  forces  hghting  against  insurgent  and  terrorist 
groups.  The  US  needs  to  adapt  to  this  type  of  warfare  by  developing  methods  and 
weapons  to  effectively  eliminate  enemy  targets  without  injuring  or  killing  innocent 
bystanders. 

1.1.2  Munitions 

Precision-guided  munitions  have  allowed  the  US  military  to  strike  specihc  targets 
with  great  accuracy.  These  bombs  and  missiles  are  dropped  from  high  altitudes  and 
are  guided  within  a  few  feet  of  the  target.  Although  such  munitions  reduce  collateral 
damage  by  reducing  the  miss  distance  from  the  actual  target,  metal  fragments  from 
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the  munition  case  travel  large  distances  and  kill  innocent  civilians  near  the  detonation 
location. 

Most  conventional  munitions  are  composed  of  a  metal  case  that  encloses  the  ex¬ 
plosive  material  and  other  components.  The  purpose  for  the  metal  shell  is  to  allow 
the  munitions  to  penetrate  tanks,  buildings,  bunkers,  and  other  material  without  dis¬ 
integrating  before  detonation.  When  the  munition  detonates,  the  internal  pressure 
creates  tensile  stress  on  the  case  as  it  expands  causing  it  to  crack  and  creating  many 
metal  fragments.  “The  sizes  of  fragments  are  a  function  of  the  rate  of  cylinder  expan¬ 
sion”  due  to  the  explosion,  which  can  be  calculated  by  Mott’s  equations  [4].  Chapter 
27  of  Cooper’s  textbook  [4]  provides  equations  that  calculate  the  velocity  and  position 
of  the  fragments.  It  can  be  calculated  that  fragments  can  travel  at  speeds  close  to 
Mach  3  and  travel  hundreds  of  feet  depending  on  the  size  of  the  fragment. 

The  need  for  munitions  with  no  to  little  fragmentation  is  growing  due  to  current 
and  forthcoming  wars.  For  that  reason,  the  Air  Force  Research  Laboratory  (AFRL) 
Munitions  Directorate  is  exploring  and  testing  munitions  that  do  not  produce  metal 
fragmentation  during  detonation.  Reducing  fragmentation  in  munitions  can  be  ac¬ 
complished  by  using  a  composite  case  made  from  materials  like  wound-carbon-fiber. 
During  the  explosion  the  high  temperatures  and  pressures  disintegrate  the  carbon- 
hber  casing,  thus  not  producing  any  fragments.  These  carbon-hber  casing  warheads 
are  a  solution  to  reducing  fragmentation.  The  next  step  in  helding  such  munitions 
is  to  develop  the  best  explosive  formulation  that  would  maximize  the  performance  of 
new  carbon-hber  casing  munitions. 
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1.2  Research  Problem  Description 


1.2.1  Introduction 

Dr.  Crochet  et  al  have  performed  numerical  analysis  on  multi-phase  explosives 
simulations  composed  of  metal  particulates  and  explosive  compounds  [7].  These  nu¬ 
merical  analyses  have  provided  utile  information  about  multi-phase  explosives.  The 
addition  of  metal  particulates  in  explosive  mixtures  increases  density  to  the  shock 
wave,  causing  a  higher  pressure  in  the  shock.  These  high-pressure  shocks  are  devas¬ 
tating  and  will  incapacitate  adversaries  in  physical  proximity  to  the  detonation  of  the 
munitions.  The  combination  of  low  fragmentary  munition  casings  and  the  increase 
in  shock  wave  pressure,  through  metalized  explosives,  has  the  potential  to  effectively 
eliminate  threats  in  modern  asymmetrical  combat  theaters. 

1.2.2  Background  and  Issues  of  Metalized  Explosives 

The  issue  with  such  a  weapon  is  that  the  physics  of  multi-phase  explosives  are 
not  well  understood.  The  inclusion  of  metal  particulates  in  explosives  represents 
a  particularly  computationally  challenging  physical  environment,  in  that  such  flows 
are  “two-phase”  flows.  In  “two-phase”  flows,  gas  physics  and  solid  particle  physics 
must  be  simultaneously  simulated.  Since  the  two  phases  represent  different  states 
of  matter,  matching  of  physical  conditions  for  the  two  different  phases  is  difficult  to 
achieve  computationally. 

In  a  simple  explosive  mixture  there  are  two  components,  the  solid  explosive  grains 
and  the  reaction  gas  products  [5].  Piston-impact  simulations  are  commonly  used  to 
determine  properties  of  the  explosive  being  tested.  Detonation  velocity  (D),  pres¬ 
sure  (P),  density  (p),  temperature  (T),  and  material  sensitivity  are  a  few  of  the 
properties  resulting  from  piston-impact  simulations.  Other  important  properties  that 
provide  an  overview  of  the  explosive  mixture  is  Run2Det  distance  and  time.  The  de- 
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scription  “Run-to-detonation”  is  abbreviated  to  “Run2Det”  in  this  report.  Run2Det 
time  is  the  point  when  the  chemical  reaction  of  the  detonation  reaches  steady-state 
and  its  detonation  properties  remain  constant.  Run2Det  distance  is  the  distance  from 
the  piston  at  which  the  shock  wave  reaches  steady-state.  As  the  impact-induced  det¬ 
onation  is  initiated,  a  shock  front  compresses  and  heats  the  explosive,  which  initiates 
an  instant  exothermic  chemical  reaction.  The  energy  released  from  the  exothermic 
reaction  continues  to  feed  the  shock  and  drives  the  shock  forward  [4].  Conventional 
explosives  composed  of  a  simple  mixture  are  computationally  simpler  because  the  ex¬ 
plosive  transitions  from  a  solid  to  a  gas  phase.  Simple  mixture  explosives  require  fewer 
equations  that  model  the  conservation  of  mass,  energy,  and  momentum  relationships. 

Simulations  become  more  complex  and  expensive  with  explosives  that  contain 
multiple  solid  components.  Multi-component  explosives  are  a  three-phase  system; 
one  phase  consists  of  a  mixture  of  condensed  explosive  and  metal  grains,  another 
phase  consists  of  the  reaction  zone,  and  the  product  phase  consists  of  solid  metal 
oxide  and  vaporized  metal  [5].  Baer  and  Nunziato  (BN)  formulated  a  system  of 
equations  that  allow  basic  simulations  of  “two-phase”  flows,  more  details  of  which 
can  be  found  in  in  Chapter  III  of  this  thesis  [3]  .  This  system  of  equations  considers 
separate  mass,  momentum,  and  energy  balance  equations  for  each  component  and 
additional  expressions  related  to  the  evolution  of  the  solid  volume  fractions  [5]. 

1.2.3  Research  Assumptions  and  Limitations 

The  complexity  of  computationally  modeling  explosive  reactions  pushes  researchers 
to  make  many  assumptions  in  order  to  obtain  approximations  of  the  reaction.  With¬ 
out  assumptions  these  simulations  would  be  almost  impossible  to  model  correctly. 
Most  of  the  modeling  methods  use  Euler’s  conservation  equations,  which  assume  that 
the  flow  is  inviscid  and  adiabatic.  Chapter  III  explains  the  assumptions  that  were 
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made  in  order  to  obtain  appropriate  equations  of  state  for  the  explosive.  Some  of 
those  assumptions  are  that  the  model  will  conserve  mass,  momentum,  and  energy, 
and  that  the  flow  is  inviscid  and  adiabatic. 

Although  a  Three-Dimension  (3-D)  simulation  would  provide  results  closer  to  ex¬ 
perimental  explosions,  it  would  be  more  computationally  and  analytically  intensive. 
Most  explosive  modeling  is  limited  to  One-Dimension  (1-D)  because  it  is  less  com¬ 
putationally  expensive  and  will  provide  a  simple  wave  solution  that  is  sufficient  to 
determine  detonation  properties. 

Certain  conditions  in  multi-phase  explosives  such  as  the  transition  from  solid  state 
to  gas  state  of  metal  particulates  can  complicate  computational  simulations.  Model¬ 
ing  methods  for  heterogeneous  explosives  can  become  very  complicated.  The  choice 
of  technique  is  dependent  upon  the  required  accuracy  and  computational  resources 
available.  Chapter  III  includes  a  detailed  description  of  the  multi-phase  modeling 
method  chosen  for  this  research. 

1.3  Research  Questions 

The  main  focus  of  this  research  is  the  study  of  the  effects  of  aluminum  in  high  ex¬ 
plosives  such  as  PBX9501  using  the  hydrocode  developed  by  AFRL.  A  hydrocode  is 
a  computational  program  for  modeling  the  behavior  of  fluid  flow  like  explosive  flows. 
This  research  utilizes  the  MPEXS  hydrocode  which  was  developed  by  Dr.  Crochet, 
a  computational  researcher  at  the  Munitions  Directorate  of  AFRL.  Dr.  Crochet’s 
numerical  scheme  demonstrated  that  detonation  velocity  has  a  strong  dependence  on 
various  parameters  such  as  the  metal  grain  size  and  metal  mass  fraction  that  is  mixed 
with  high  explosive  HMX  {C^HsN^Os)  [5].  Al  grain  size  and  initial  metal  mass  frac¬ 
tion  (Am)  will  be  the  independent  variables  of  interest  to  study  the  detonation  effects. 
The  MPEXS  code  will  be  used  to  evaluate  the  effects  of  aluminum  particulates  in  the 
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high  explosive  PBX9501  and  provide  a  recommendation  on  an  explosive  mixture. 

Determining  the  effects  of  Al  grain  size  and  mass  fraction  on  detonation  properties 
for  PBX9501  are  the  main  focus,  but  the  MPEXS  hydrocode  will  also  be  evaluated. 
The  MPEXS  code  will  be  evaluated  by  performing  several  studies  to  determine  the 
stability  and  limitations  of  the  code.  The  MPEXS  code  needs  to  output  comparable 
data  to  other  hydrocodes.  The  stability  of  the  MPEXS  code  will  be  determined  by 
looking  for  a  convergence  in  the  data  as  the  grid  resolution  is  increased.  Performing 
a  wide  range  of  simulations  will  insure  that  minimums  and  maximums  will  be  tested, 
providing  an  idea  of  the  limitations  of  the  code. 

1.4  Overview  of  Subsequent  Sections 

Chapter  I  includes  a  background  on  current  asymmetrical  warfare  and  the  prob¬ 
lems  presented  to  the  US  and  its  allies.  One  of  the  problems  faced  is  the  high  proba¬ 
bility  of  injury  or  death  to  innocent  civilians  resulting  from  conventional  munitions. 
A  solution  to  such  problem  is  the  development  of  low  fragmentary  munition  casings 
and  explosive  formulations  that  yield  strong  blast  effects  within  a  small  radius  and 
then  rapidly  drop  off  in  strength.  Chapter  I  also  includes  the  research  problems  and 
the  research  questions  to  be  answered.  Chapter  II  will  present  a  literature  review  on 
the  fundamentals  of  detonation  theory,  previous  research  on  the  effects  of  aluminum 
on  detonation  velocity,  and  multi-phase  explosive  modeling  techniques.  Chapter  III 
will  present  a  detailed  methodology  of  the  simulations  and  data  collection  procedures 
that  was  implemented  to  obtain  the  sought  after  results.  Chapter  IV  will  present  the 
results  of  the  simulations  that  were  performed.  Finally,  Chapter  V  will  summarize 
the  data,  and  make  conclusions  about  the  effects  of  aluminum  on  PBX9501  and  the 
performance  of  the  MPEXS  code. 


II.  Literature  Review 


2.1  Introduction 

To  understand  the  effects  metal  particulates  have  in  explosive  materials,  one  must 
begin  with  the  basics.  This  chapter  will  include  a  literature  review  of  the  fundamentals 
of  detonation  to  understand  the  physics,  thermodynamics,  and  chemical  reactions 
that  occur  in  an  explosive  detonation.  A  shock  wave  traveling  radially  out  from  the 
center  of  the  explosion  provides  very  useful  information  of  the  detonation.  Detonation 
velocity,  pressure,  and  temperature  are  a  few  of  the  parameters  that  are  studied. 
Throughout  the  years  much  of  empirical  data  has  been  collected  on  explosives  and 
equations  have  been  developed  based  on  this  data.  In  recent  years  new  heterogeneous 
explosives  have  emerged.  These  explosives  might  be  a  mixture  of  several  explosive 
materials  or  a  mixture  of  explosive  materials  and  metalized  powders.  This  literature 
review  also  includes  current  simulation  work  on  heterogeneous  explosive  materials,  as 
well  as  explosive  materials  with  aluminum  particulates  which  is  the  area  of  interest 
for  this  thesis. 

2.2  Fundamentals  of  Detonation  Theory 

2.2.1  Simple  Theory 

An  explosive  detonation  is  a  very  complicated  fluids  problem  because  it  is  a  shock 
wave  traveling  out  radially  with  a  rapid  exothermic  chemical  reaction  zone  occurring 
behind  the  shock  front  [4] .  By  considering  a  simple  model  of  detonation,  it  puts  this 
complicated  fluids  problem  into  simpler  terms.  The  “simple  theory”  or  “Zeldovich, 
Von  Neumann,  and  Deering  (ZND)  theory”  mathematically  simplifies  detonation  by 
quantifying  thermodynamic  and  kinetic  terms  into  first-order  engineering  problems 
[4].  The  ZND  theory  makes  the  following  six  assumptions  [4]: 
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First:  The  flow  is  1-D  (laminar). 


Second:  The  detonation  front  is  a  discontinuity  because  there  is  an  instantaneous 
jump  in  the  pressure. 

Third:  All  the  reaction  gases  are  in  chemical  and  thermodynamic  equilibrium. 

Fourth:  The  chemical  reaction  zone  is  infinitely  thin. 

Fifth:  The  detonation  velocity  is  constant  and  therefore  the  products  leave  at 
the  same  state  and  are  independent  of  time. 

Sixth:  The  detonation  products  can  change  after  leaving  the  reaction  zone. 

Based  on  these  assumptions,  a  detonation  is  defined  as  a  shock  wave  passing 
through  the  explosive  (explosive  or  energetic  is  interchangeably  through  out  this  the¬ 
sis)  with  a  following  rarefaction  wave  or  gas  expansion.  This  wave  is  known  as  a 
Taylor  Wave,  named  after  the  British  mathematician  Geoffrey  Ingram  Taylor  who 
developed  the  Equation  of  State  (EOS)s  for  the  wave.  The  Taylor  wave  can  be  seen 
in  Figure  1.  “The  shock  front,  chemical  reaction,  and  the  leading  edge  of  the  rarefac¬ 
tion  are  all  in  equilibrium;  so  they  are  all  moving  at  the  same  speed,  which  we  call 
detonation  velocity  (D)  [4].”  The  steady-state  detonation  condition  for  the  products 
is  known  as  the  Chapman-Jouguet  (CJ)  point,  shown  by  the  blue  dot  in  Figure  1. 
The  Von  Neuman  Spike  is  the  necessary  energy  to  activate  the  chemical  reaction, 
see  Figure  1.  The  distance  between  the  vertical  dashed  line  and  vertical  solid  line  in 
Figure  1  represents  the  reaction  zone  of  the  detonation  wave.  For  ZND  theory,  the 
reaction  zone  is  assumed  to  be  infinitely  thin  or  zero  as  stated  by  the  forth  assumption 
listed  above.  Another  important  detonation  parameter  is  the  CJ  pressure  {Pcj)- 
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2.2.2  Wave  Properties 


In  explosives,  piston-impact  simulations  are  commonly  utilized  to  determine  prop¬ 
erties  of  the  explosive  materials  being  tested.  A  piston-impact  problem  is  simulated 
by  a  rigid  piston  moving  with  constant  Up  inside  a  rigid  tube  which  contains  explo¬ 
sive  materials  and  gases.  The  denotation  is  initiated  when  the  piston  compresses  the 
explosive  materials  and  gases.  A  forward-moving  shock  wave  travels  with  velocity 
D.  To  understand  this  piston-impact  problem,  the  flow  held  between  the  piston  and 
the  shock  must  be  solved.  D  is  an  important  parameter  that  could  be  solved  by 
using  the  Hugoniot  relationships.  A  suitable  how  connecting  the  piston  and  the  hnal 
state  must  be  found  [9].  The  solution  can  consist  of  a  rarefaction  wave.  When  the 
piston  velocity  is  less  than  detonation  velocity,  the  front  of  the  wave  still  exists  and 
a  rarefaction  wave  is  needed  “to  reduce  the  velocity  at  the  front  to  that  at  the  piston 
[9].”  Figure  1  shows  the  detonation  wave  front  and  the  rarefaction  wave  that  follows. 
Figure  1  is  based  on  a  hgure  by  Fickett  and  Davis  [9]. 


Figure  1.  Detonation  Wave 
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The  final  state  as  described  by  Fickett  and  Davis,  is  where  the  flow  stabilizes 
just  after  the  reaction  zone  of  the  detonation  wave  [9].  “The  conservation  conditions 
require  that  the  final  state  point  in  the  pressure  vs  specific  value  plane  lie  on  both  this 
Hugoniot  curve  and  the  Rayleigh  line  [9].”  The  final  state  is  given  by  the  intersection 
of  the  Rayleigh  line  and  the  Hugoniot  curve  as  shown  in  Figure  2.  Figure  2  is  based  on 
a  figure  by  Fickett  and  Davis  [9].  There  are  three  different  final  state  solutions.  1)  If 
D  is  less  than  the  CJ  detonation  velocity  (Dcj),  the  red  Rayleigh  line  does  not  cross 
the  blue  Hugoniot  curve,  therefore  there  is  no  steady-state  solution  [9].  2)  If  D  is  equal 
to  Dqj,  the  red  Rayleigh  line  crosses  the  blue  Hugoniot  curve  at  one  location,  the 
CJ  point,  therefore  there  is  one  steady-state  solution  [9].  3)  If  D  is  greater  than  the 
Dcj,  the  red  Rayleigh  line  crosses  the  blue  Hugoniot  curve  at  two  locations,  therefore 
there  are  two  possible  solutions  [9].  The  upper  or  strong  intersection  where  the  flow 
is  subsonic  is  utilized  for  steady  state  analysis  [9].  The  lower  or  weak  intersection  is 
rejected  because  the  flow  is  supersonic  [9]. 


Hugoniot  Curve 


Figure  2.  Hugoniot  Curve  and  Rayleigh  Lines  Relationships:  Final  State  Conditions 
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2.2.3  Solving  for  Explosive  Parameters 


There  are  a  few  important  parameters  that  need  to  be  determined  for  explosive 
analysis  and  those  are  CJ  density  (pcj),  CJ  particle  velocity  {ucj),  Pcj,  and  D. 
Empirical  data  shows  that  there  is  a  relationship  between  initial  density  (po)  and  the 
Pcj-  When  plotting  the  empirical  data  on  a  log  scale  it  is  seen  that  there  is  a  fairly 
linear  relationship.  Therefore,  Cooper  derived  Equation  1  to  solve  for  pcj  [4]. 

PCJ  =  l.SSQp^-^^  [4]  (1) 

where  pcj  is  CJ  density,  and  po  is  initial  density. 

The  Pcj  and  ucj  can  be  determined  from  initial  density  and  D  empirical  values 
using  mass  and  momentum  conservation  equations.  Equation  2  provides  an  estimate 
of  Pcj  within  5%  of  experimentally  measured  values  [4].  Equation  3  demonstrates 
how  Ucj  can  be  calculated  when  Pcj  and  po  are  known  [4]. 

Pcj  =  PoD\l  -  0.7125pr')  [4]  (2) 


Ucj 


Pcj 

PoD 


[4] 


(3) 


where  Pcj  is  CJ  pressure,  po  is  initial  density,  D  is  detonation  velocity,  and  ucj  is 
CJ  particle  velocity. 

The  Equations  2  and  3  provide  acceptable  approximations  for  certain  explosive 
parameters  at  the  CJ  state.  It  is  also  important  to  determine  explosive  parameters  at 
other  conditions  along  the  Hugoniot  curves.  Hugoniots  are  empirically  derived  EOS 
relating  six  different  variable  pairs  such  as  pressure  vs  shock  velocity  {P  —  U),  pres¬ 
sure  vs  particle  velocity  {P  —  u),  pressure  vs  speciffc  volume  {P  —  v),  shock  velocity 
vs  particle  velocity  {U  —  u),  shock  velocity  vs  speciffc  volume  {U  —  v),  and  parti- 
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cle  velocity  vs  specific  volume  {u  —  v)  [4],  Estimating  Hugoniot  relationships  can 
be  approximated  more  accurately  using  computer  codes  that  use  empirically  derived 
nonlinear  EOS.  Due  to  the  complexity  of  the  nonlinear  EOS  these  codes  can  be  com¬ 
putationally  expensive  requiring  large  computers  or  codes  are  not  readily  available 
to  most  engineers.  However,  there  are  simple  empirical  correlations  that  will  provide 
solutions  with  reasonable  accuracy  based  on  P/Pcj  ratios.  The  data  correlates  to 
two  regions;  for  reduced  pressure  ratios  above  0.08  use  Equation  4  and  for  reduced 
pressure  ratios  below  0.08  use  Equation  5  [4]. 


p 

Pcj 


>  0.08: 


P 


2.142  -  1.7315 


(4) 


p 

Pcj 


<  0.08: 


P 


=  235 


u 


UCJ 


-8.71 


[4] 


(5) 


where  P  is  pressure,  Pcj  is  CJ  pressure,  ucj  is  CJ  particle  velocity,  and  particle 
velocity  (u). 

There  are  several  methods  to  develop  equations  of  state  for  an  explosive  mate¬ 
rial.  Equations  of  state  for  explosive  products  are  only  engineering  approximations 
for  simulating  explosive  behaviors  [9].  As  discussed  above,  one  of  the  methods  is  to 
End  a  Hugoniot  relationship  between  the  six  different  variable  pairs.  These  relation¬ 
ships  can  be  done  with  nonlinear  EOS  using  computers  or  using  simple  empirical 
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correlations.  These  simple  theory  EOSs  have  many  limitations  that  are  based  from 
the  assumptions  made  in  the  Simple  Theory  section.  So  far  we  have  focused  on  the 
EOS  for  the  explosive  material,  but  now  we  need  to  look  for  EOS  for  the  reaction 
products.  Fickett  and  Davis  discussed  several  forms  that  have  been  utilized  to  deter¬ 
mine  reaction  product  EOS  [9].  Some  of  these  forms  account  for  the  chemistry  in  the 
reaction  and  some  that  do  not.  There  are  four  EOS  without  chemistry  discussed  by 
Fickett  and  Davis;  1)  the  constant-7  form,  2)  the  constant-/?  equation  of  state,  3)  the 
constant-a  equation  of  state,  and  4)  the  Jones-Wilkins-Lee  (JWL)  equation  of  state 
[9] .  One  of  the  simplest  form  is  the  constant-7  form  which  yields  Equation  6  for  Pcj 
and  Equation  7  for  Dcj. 


Pcj  =  2(7  -  l)p,q  [9]  (6) 

Dcj^  =  2(7'  -  l)q  [9]  (7) 

where  Pcj  is  CJ  pressure,  Dcj  is  detonation  velocity,  heat  capacity  ratio  (7)  =  3,  po 
is  initial  density,  and  specihc  heat  (g). 

Another  form  is  the  constant-/?  equation  of  state  were  /?  is  constant  for  all  pressures 
and  specihc  volumes.  The  equations  for  Pcj  and  Dcj  are  complicated,  therefore 
Fickett  and  Davis  derived  Equations  8  and  9  that  approximate  Pcj  and  Dcj  based 
on  empirical  correlations  for  the  explosive  Composition  B  [9]. 

Pcj  cx  [9]  (8) 

Dcj^  oc  pP\P^  [9]  (9) 

where  Pcj  is  CJ  pressure,  Dcj  is  detonation  velocity,  po  is  initial  density,  and  q  is 
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specific  heat. 

Another  form  is  the  constant-a  equation  of  state.  These  EOS  are  similar  to  the 
constant-/?  equation  of  state,  but  a  is  a  function  of  pressure,  where  as  /?  is  a  function 
of  specihc  volume.  Just  like  the  constant-/?  equation  of  state,  the  Pcj  and  Dqj 
equations  are  complicated  and  require  a  computer  code  to  solve.  Another  form  of  the 
constant-/?  EOS  is  the  JWL  equation  of  state,  but  it  uses  the  Gruneisen  coefficient 
instead  of  /?.  The  four  different  forms  can  be  utilized,  but  according  to  Eickett  and 
Wood  the  constant-/?  equation  of  state  has  the  most  reasonable  properties  in  the 
region  of  interest  [9]. 

Once  a  form  of  the  EOS  has  been  chosen,  important  parameters  like  P  and  D 
can  be  determined.  The  calculated  P  and  D  values,  which  were  derived  from  the 
EOS  must  be  compared  to  empirical  data.  This  comparison  will  determine  if  the 
coefficients  utilized  for  the  P  and  D  equations  provided  a  good  correlation  to  the 
empirical  data.  If  the  coefficients  did  not  provide  a  good  correlation  then  they  must 
be  reconsidered.  The  EOS  should  be  calibrated  using  the  empirical  data  until  a 
stronger  correlation  is  derived. 

2.3  Heterogeneous  Explosive  Materials 

As  previously  stated,  determining  explosive  parameters  of  pure  explosive  can  be 
accomplished  with  certain  conhdence,  but  it  can  get  computationally  intensive  de¬ 
pending  on  the  accuracy  needed  and  the  modeling  methods  utilized.  This  problem 
is  exacerbated  when  the  explosive  materials  are  heterogeneous.  Heterogeneous  im¬ 
plies  the  explosive  mixture  is  composed  of  different  materials,  phases,  or  grains  of 
different  densities.  It  is  important  that  heterogeneous  explosive  materials  are  studied 
because  almost  all  explosives  utilized  in  the  military  or  commercially  are  heteroge¬ 
neous.  Some  explosives  use  binding  materials  to  hold  the  explosive  crystals  together 
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as  well  as  to  decreases  the  shock  sensitivity  of  the  explosive  by  preventing  the  ex¬ 
plosive  grains  from  touching  each  other,  making  it  safer  to  handle.  Polymer-Bonded 
Explosive  (PBX)  and  Livermore’s  High-Energy  Explosive  (LX)  are  examples  of  ex¬ 
plosives  that  use  “plastic”  as  a  binder  material.  Other  heterogeneous  explosives  mix 
two  or  more  explosive  materials,  for  example  Composition  B  consists  of  a  mixture  of 
Research  Department  Explosive  (RDX)  and  Trinitrotoluene  (TNT). 

The  microstructure  of  these  heterogeneous  materials  introduce  internal  boundary 
effects  that  change  the  thermal/mechanical/chemical  behavior  of  the  explosive  [2]. 
The  microscale,  at  the  atomic  level,  is  usually  not  utilized  for  these  types  of  exper¬ 
iments  because  it  requires  greater  computational  resources.  To  get  an  idea  between 
the  micro  and  meso  scales.  Dr.  Baer  gives  an  example  that  a  1pm  long  solid  material 
flake  would  contain  10^^  atoms  [2].  Therefore,  the  mesoscale  is  more  appropriate  for 
these  explosive  experiments. 

Some  of  Dr.  Baer’s  prior  work  has  demonstrated  that  particle  diameter  and  mor¬ 
phology  are  very  important  parameters  to  the  formulation  of  explosive  materials  [2]. 
Many  performance  parameters  like  shock  sensitivity,  pressure,  detonation  velocity, 
and  particle  velocity  are  believed  to  be  directly  affected  by  particle  size  and  mor¬ 
phology  of  the  explosive.  In  order  to  trigger  a  chemical  reaction  a  certain  pressure 
threshold  must  be  met.  Tests  indicated  that  hue  grain  HMX,  with  average  particle 
diameter  of  10-15pm,  required  higher  pressures  to  trigger  the  reaction  than  other 
larger  grain  HMX  [2].  Most  of  the  current  work  is  either  done  through  computer 
simulation  or  by  impact  loading  experiments,  but  better  “experimental  techniques 
for  measurement  at  the  mesoscale  crystal  level  are  currently  in  development  [2].” 
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2.3.1  Uniform  Mixing  Methods 


Heterogeneous  materials  when  mixed  properly  should  have  uniform  mixture  of 
all  the  particles.  The  explosive  mixture  should  be  uniform  enough  that  samples  of 
the  same  mixture  should  not  vary  from  sample  to  sample.  As  with  real  explosives, 
computer  simulations  of  explosives  should  have  uniform  mixture.  The  heterogeneous 
material  must  be  modeled  by  using  a  computerized  method  that  randomly  distributes 
the  varied  geometry  crystals  into  a  closely  packed  configuration,  so  that  the  mixture 
is  as  uniform  as  possible  [2].  Due  to  the  infinite  number  of  possible  combinations,  the 
size  and  geometry  of  the  explosive  crystals  must  be  limited  to  a  certain  range  of  sizes 
and  geometries.  Since  the  main  concern  is  particle/crystal  size,  the  geometry  can 
be  limited  to  one  shape,  even  though  in  real  explosive  materials  a  myriad  of  crystal 
shapes  exist. 

Dr.  Baer  mentions  two  computer  simulation  methods  to  randomly  distribute  and 
pack  particles  of  a  heterogeneous  material  [2].  One  of  those  methods  is  the  drop- 
and-fill  method  were  randomly  sized  particles  are  dropped  into  a  3-D  space,  then  a 
simulated  gravitational  force  will  settle  the  particles  at  the  bottom  of  the  3-D  space. 
The  problem  with  this  method  is  that  “the  infiuence  of  “gravity” ,  leads  to  interbody 
collision  and  contact  frictional  effects  [2].”  Another  method  to  randomly  distribute 
and  pack  particles  of  arbitrary  sizes  is  by  using  a  combination  of  Monte  Carlo  (MC) 
and  Molecular  Dynamics  (MD)  methods.  The  MC  method  uses  a  random  number 
generator  based  on  a  chosen  probability  distribution  function  to  determine  different 
particle  sizes.  Figure  3  is  based  on  a  ffgure  by  Dr.  Baer  [2].  The  randomly  sized 
particles  from  the  MC  method  are  placed  in  a  3-D  space  and  distributed  randomly, 
as  shown  in  Figure  3a).  An  MD  simulation  method  places  particles  by  giving  them  a 
velocity  and  and  allowing  them  to  collide  as  shown  in  Figure  3b).  This  collision  dy¬ 
namics  ensure  that  the  particles  mix  with  each  other,  creating  a  random  distribution. 
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In  order  to  pack  them  to  a  desired  volume  fraction,  the  density  inside  the  3-D  space 
is  increased  dynamically  until  the  particles  are  locked  in  place,  Figure  3c)  [2], 


a)  Distribute  Randomly  b)  Collision  Dynamics  c)  Compress  Space 
Figure  3.  Monte  Carlo  Method  for  Randomly  Packing  Particles 

2.3.2  Experimental  Measuring  Techniques 

As  mentioned  earlier,  computational  simulations  are  effective  and  provide  useful 
approximation,  but  experimental  testing  provide  more  realistic  results.  Experimen¬ 
tal  techniques  for  measurement  at  the  mesoscale  are  currently  in  development  by 
organization  such  as  Sandia  National  Laboratories  (SNL).  A  new  technique  being 
explored  at  SNL  is  a  line-imaging  Optically  Recording  Velocity  Interferometry  Sys¬ 
tem  (ORVIS)  to  provide  data  on  particle  velocity  in  shock-compressed  heterogeneous 
materials  [2],  ORVIS  uses  electromagnetic  waves  to  measure  properties  in  the  shock 
wave  and  produces  a  3-D  line-image  of  the  particle  velocity  versus  time  and  distance. 
The  time  resolution  of  the  ORVIS  technique  is  approximately  3  nanoseconds  and  a 
spatial  resolution  of  approximately  3  micrometers.  For  examples  of  an  ORVIS  particle 
velocity  measurement  of  a  compaction  wave  see  work  by  Dr.  Baer  [2]. 

Figure  4a  shows  a  schematic  of  the  line-imaging  ORVIS  based  on  a  hgure  from  Dr. 
Ao  [1].  The  light  beam  from  a  moving  object,  such  as  a  shock  wave,  is  reflected  to  a 
series  of  mirrors  and  optics  to  collimate  the  light.  The  light  signal  then  pases  through 
a  splitter,  where  the  signal  either  goes  to  mirror  one  (Ml)  or  mirror  two  (M2).  The 
M2  signal  gets  delayed  by  a  delay  etalon  and  the  Ml  signal  continues  with  no  delay. 
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Both  the  signals  meet  at  a  crossing  plane,  shown  in  Figure  4a  and  b.  The  interfering 
fringes  get  translated  into  an  output  signal  of  velocity  and  time  [1]. 


light  beam 
from  object 


(b)  Crossing  Plane 


Figure  4.  ORVIS  Schematic  and  the  Crossing  Plane 


The  ORVIS  experimental  technique  opens  many  possibilities  in  exploring  hetero¬ 
geneous  materials,  but  there  are  a  few  issues  with  the  technique  [2].  One  of  the  issues 
is  the  large  amount  of  data  that  is  collected  needs  to  be  stored,  analyzed,  and  inter¬ 
preted.  SNL  is  looking  at  ways  to  improve  the  post  processing  of  large  amounts  of 
data.  Despite  the  minor  drawbacks,  this  technique  shows  incredible  promise  in  allow¬ 
ing  for  the  capture  of  useful  data  from  physical  tests,  providing  a  practical  alternative 
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to  reliance  on  compntational  techniqnes.  The  combination  of  simnlation  and  experi¬ 
mental  data  will  allow  for  a  better  nnderstanding  of  the  thermal/mechanical/chemical 
behaviors  of  heterogeneons  explosives. 

2.3.3  Explosive  Materials  with  Aluminum  Powders 
2.3.3. 1  Inert  Aluminum 

The  addition  of  metal  particulates  in  explosives  is  not  a  new  concept,  but  there 
remains  a  lack  of  understanding  on  what  occurs  in  the  reaction.  Several  researchers 
have  developed  EOS  for  metalized  explosives  in  an  attempt  to  computationally  model 
the  behavior.  Dr.  Crochet  developed  a  simulation  modeling  code,  MPEXS,  to  study 
the  effects  of  metalized  explosives  [5].  In  particular,  he  looked  at  the  high  explosive 
HMX  and  Al  by  varying  initial  metal  mass  fraction  (A^)  by  15%,  25%,  50%  for  Al 
grain  mean  diameter  {dm)  of  150/im  and  lOOnm.  As  A^  increases,  the  amount  of  solid 
explosive  available  for  combustion  decreases,  therefore  affecting  detonation  properties 
like  detonation  velocity,  pressure,  and  temperature.  For  inert  metalized  explosives, 
as  Xm  increases  the  heat  released  from  the  compressive  work  behind  the  detonation 
wave  decreases  [5].  The  heat  transfer  from  the  gas  phase  is  mostly  influenced  by  the 
combustion  reaction  energy,  which  is  dependent  on  the  aluminum  particle  size,  i.e. 

dm  [5] . 

The  results  for  detonation  wave  speed,  end-state  pressure,  and  end-state  temper¬ 
ature  for  the  three  different  metal  mass  fractions  at  various  metal  particle  diameters 
are  presented  in  Dr.  Crochet  dissertation  [5].  It  was  concluded  by  Dr.  Crochet  that 
the  detonation  wave  speed  increases  as  the  mass  fraction  decreases  [5].  Therefore,  the 
amount  of  Al  present  in  the  mixture  directly  affects  the  wave  speed.  As  the  mass  frac¬ 
tion  decreases,  both  pressure  and  temperature  increase.  The  wave  speed,  pressure, 
and  temperature  remain  constant  for  small  aluminum  grain  sizes  until  a  threshold  is 
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met,  then  it  increases.  Dr.  Crochet  suggests  that  “the  influence  of  size-dependent  re¬ 
laxation  processes  approach  a  threshold,  beyond  which  more  rapid  phase  equilibrium 
has  a  negligible  effect”  on  wave  speed,  pressure,  and  temperature  [5]. 

Dr.  Crochet’s  inert  metal  model,  when  compared  to  the  experimental  data  of 
Gogulya  et  al.  [11],  under  predicted  the  detonation  wave  velocities  for  aluminum 
grain  sizes  of  20,  50,  and  150  microns  [5].  One  of  the  factors  that  is  believed  to 
reduce  the  detonation  wave  velocity  is  that  as  the  \°^  increases,  the  explosive  gets 
diluted  and  combustion  energy  available  decreases  [5].  The  other  factor  is  as  the 

increases,  gas-phase  energy  decreases  [5].  These  decreases  in  energy  reduce  the 
detonation  wave  velocities.  For  the  500  nm  aluminum  grain  size,  the  inert  metal 
model  over  predicted  the  detonation  wave  velocity,  “possibly  due  to  the  simplified 
form  of  the  drag  and  heat  transfer  coefficients”  [5]. 

2. 3. 3. 2  Reactive  Aluminum 

Dr.  Crochet  also  performed  an  analysis  of  reactive  metal  in  the  explosive  by 
varying  two  parameters;  the  metal  grain  size  dm  and  oxidation  rate  {Kpr)  [5].  Based 
on  the  metal  oxidations  laws  utilized  in  the  simulations,  there  were  three  values  of 
Kpr  tested.  Oxidation  rates  of  Kpr=  1.9x10“®  w? / s  and  Kpr=  1.9x10“®  rri^ / s  were 
taken  from  suggested  values  based  on  Fedorov  and  Kharlamova’s  research  [5].  The 
other  oxidation  rate  is  that  of  inert  aluminum,  Kpj.=  0  m? / s  [5].  The  results  of  the 
detonation  wave  velocity  with  the  two  varying  parameters,  dm  and  Kpr  are  presented 
in  Dr.  Crochet  dissertation  [5].  For  large  aluminum  particles,  dm  greater  than  5 
microns,  the  oxidation  prefactor  seems  not  to  affect  the  detonation  velocity.  Due  to 
a  more  rapid  metal  oxidation  and  heat  transfer,  it  is  observed  that  the  detonation 
velocity  increases  as  Kp^  increases  for  dm  less  than  5  microns  [5].  Similar  to  the  inert 
metal,  there  is  a  threshold  of  metal  particle  size,  above/below  which  the  size  does  not 


22 


affect  the  detonation  velocity. 

The  study  performed  by  Dr.  Crochet  consisted  on  the  explosive  HMX  and  alu¬ 
minum,  where  aluminum  was  either  inert  or  reactive  [5].  In  both  of  these  cases,  it  was 
observed  that  aluminum  particle  size  had  a  direct  impact  on  many  explosive  prop¬ 
erties  like  detonation  wave  velocity,  end-state  pressure,  and  end-state  temperature. 
Future  work  in  this  area  includes  improving  the  rehnement  of  spatial  grid  utilized  near 
the  detonation  wave  to  properly  model  gas  physics,  thermal  and  chemical  behaviors. 

2.3.4  An  Analysis  of  Shock-Induced  Reactions  in  FeaOs+A/  +  Teflon® 
Powder  Mixtures:  Reactive  Material  Section 

Yang  et  al  [15]  conducted  a  series  of  experiments  to  perform  a  preliminary  anal¬ 
ysis  of  shock-induced  chemical  reactions.  The  powder  mixture  for  this  series  of  ex¬ 
periments  was  a  mixture  of  aluminum,  hematite  {Fe20f),  and  Polytetrafluoroethy- 
lene  (PTFE)  also  known  as  Teflon®.  These  powder  mixture  samples  were  impacted 
by  a  metal  “flyer”,  on  impact  it  sends  a  shock  wave  through  the  powder  mixture. 
This  shock  wave  initiates  the  chemical  reaction  and  continues  to  feed  the  shock  wave. 
There  are  two  Poly vinyli dene  Fluoride  (PVDF)  pressure  gauges  located  at  different 
positions  between  the  powder  mixture  that  measure  shock  pressures  .  These  gauges 
are  utilized  to  determine  what  is  happening  in  the  chemical  reaction,  which  occurs 
under  the  high  pressure  shock  wave  loading. 

During  the  Fe203-fA/-|- Teflon®  powder  mixture  chemical  reaction,  Yang  et  al 
believe  that  two  main  reactions  occur.  One,  the  reaction  between  Fe20^  and  Al 
and  the  other  between  the  Teflon®  and  aluminum  [15].  When  aluminum  reacts  with 
Teflon®  or  hematite  it  produces  a  gas  product  of  aluminum  fluorides  (AlFj),  with  j 
representing  the  number  of  fluoride  atoms  present.  The  variation  of  shock  strength 
during  the  propagation  in  the  mixture  is  thought  to  be  determined  by  the  product 
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gas  AlFj  [15].  It  was  concluded  by  Yang  et  al  that  the  AlFj  formed  depends  on  both 
the  shock  pressure  and  temperature  in  the  mixture.  As  the  shock  pressure  increases, 
more  Fe203  and  Teflon®  react  with  aluminum;  therefore  reducing  the  amount  of 
fluoride  atoms  present  in  the  product  gas  AlFj  [15]. 

When  dealing  with  shock  induced  chemical  reactions  and  explosive  reactions,  there 
are  many  assumptions  made  due  to  the  complexity  and  speed  of  events.  In  order  to 
simplify  the  complexity  of  these  problems,  the  initial  assumptions  made  result  in 
non  desired  EOS  variables  and  coefficients.  Yang  et  al  concluded  that  the  variables 
describing  mechanical  and  thermodynamic  reactants  and  products  may  need  to  be 
improved  [15].  These  variables  change  too  fast  that  chemical  equilibriums  might  not 
be  reached,  but  it  is  difficult  to  determine  the  state  of  the  reaction.  Yang  et  al  also 
concluded  that  a  further  understanding  is  needed  on  how  to  formulate  the  conditions 
under  which  the  shock  pressure  increases  and  decreases  during  wave  propagation  [15] . 

2.4  Conclusion 

The  literature  reviewed  herein  served  as  the  foundation  for  the  research  described 
in  this  document.  The  main  goal  of  this  thesis  is  to  study  the  effects  of  aluminum 
particulates  in  certain  explosive  materials  via  a  computational  simulation.  This  chap¬ 
ter  included  a  review  of  the  fundamentals  of  detonations  of  explosives  included  the 
basic  information  of  the  important  parameters  that  are  studied  in  explosives.  Also, 
included  in  this  chapter  was  a  review  of  current  simulation  work  on  heterogeneous 
explosives  with  inert  and  reactive  aluminum  methods.  This  chapter  included  a  review 
on  computational  procedures  to  randomize  mixtures  of  particles  to  simulate  as  close 
as  possible  real  explosives. 
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III.  Methodology 


3.1  Introduction 

The  code  utilized  for  the  1-D  multi-phase  simulation  of  Al  particulates  in  explo¬ 
sives  needs  to  be  consistent  with  thermodynamic  or  entropy  conditions  that  exist  in 
the  explosive  reaction,  incorporate  the  EOS,  and  the  reactive  bnrn  models  between 
phases  dnring  the  reaction.  As  explained  earlier,  the  effects  of  metal  mass  fraction 
(Am)  and  particle  size  {dm)  are  the  parameters  of  interest.  Some  modihcations  to  the 
code  might  be  necessary  to  properly  analyze  the  effects  of  different  explosive  com- 
ponnd  formulations.  These  snggested  modihcations  will  be  implemented  in  fnture 
versions  of  the  code.  Simnlations  will  be  performed  on  several  formulations  of  explo¬ 
sives,  different  particle  grain  sizes,  and  different  metal  mass  fractions.  This  research 
will  look  at  the  stability  of  the  code  by  rnnning  several  tests  to  determine  consis¬ 
tency  and  limits  of  the  code.  This  research  will  look  at  trends  in  several  detonation 
properties  for  metalized  PBX9501.  In  this  section  a  detailed  description  of  the  mnlti- 
phase  code  will  be  presented.  Also,  the  methodology  that  will  be  ntilized  to  rnn  the 
simulations  and  how  the  data  will  be  analyzed. 

3.2  Multi-Phase  Code  Modeling 

To  properly  evaluate  the  effects  of  aluminum  metal  particulates  in  multi-phase 
explosives,  this  research  will  use  a  1-D  multi-phase  computational  numerical  method 
code,  MPEXS,  to  simulate  the  explosion  by  predicting  thermomechanical  flow  held 
properties  at  the  macroscopic  level.  The  MPEXS  code  utilized  for  this  research 
was  developed  by  Dr.  Michael  Crochet,  a  computational  researcher  at  the  Munitions 
Directorate  of  AFRL.  The  MPEXS  code  is  an  evolving  code  that  will  be  continuously 
modihed  to  provide  better  explosive  parameters  and  ease  of  use  by  Department  of 
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Defense  (DOD)  users.  The  main  goal  of  this  thesis  is  to  use  Dr.  Crochet’s  code  to 
evaluate  the  effects  of  aluminum  in  PBX9501,  a  multi-phase  explosive.  Other  goals 
are  to  test  the  accuracy  and  performance  of  the  code  and  provide  suggestions  to 
AFRL. 

This  section  will  provide  an  explanation  of  how  the  fundamental  equations  of  state 
were  derived  for  the  MPEXS  hydrocode.  It  will  provide  a  description  of  the  relevant 
input  hies  and  subroutines  hies  that  the  user  will  modify,  how  the  MPEXS  code  is 
executed,  and  the  output  hies  it  produces.  All  other  subroutines  and  equations  asso¬ 
ciated  with  the  code  are  described  in  more  detail  in  Dr.  Crochet’s  Ph.D.  dissertation 
[5].  “The  code  is  intended  to  provide  the  end  user  with  some  hexibility  in  prescrib¬ 
ing  equations  of  state,  burn  models,  and  other  constitutive  relations  by  utilizing  a 
modular  source  code  structure,  where  only  a  limited  number  of  subroutines  require 
user  interaction  [6].”  This  “beta”  version  of  the  MPEXS  hydrocode  will  be  contin¬ 
uously  modihed  until  a  hnal  version  is  user  friendly  and  provides  the  most  accurate 
simulation  results  for  metalized  explosives. 

There  are  diherent  models  that  could  be  utilized  to  simulate  explosive  reactions. 
The  One-Phase  Model  simulates  simple  explosions  and  provide  overall  how  properties 
but  no  information  is  provided  on  what  occurred  in  the  gas  or  solid  phases.  The 
Two-Phase  Model  provides  how  properties  for  both  the  solid  and  gas  phase  of  the 
explosion.  The  Two-Phase  Model  is  utilized  for  this  research  to  run  simulation  on 
simple  explosive  reactions  like  PBX9501  with  no  aluminum  particulates.  A  more 
complicated  model  like  the  Multi-Phase  Model  with  N  Solid  Components  is  utilized 
for  this  thesis  when  running  simulations  with  aluminum  particulates.  The  Multi- 
Phase  Model  allows  the  simulation  of  explosives  with  more  than  one  solid  phase,  that 
is,  either  multiple  solid  explosive  compounds  or  an  explosive  solid  compound  with  a 
solid  metal  compound.  The  following  sections  provide  a  more  detailed  description  of 
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the  different  models  to  include  the  derivations  of  the  systems  of  equations  for  each. 


3.2.1  One-Phase  Model 

The  one-phase  model  can  predict  flow  properties  of  simple  explosives  at  both  the 
macro  and  micro  scale.  Macro  scale  is  when  objects  can  be  seen  with  the  naked  eye 
and  micro  scale  is  when  objects  need  a  magnihcation  device  to  see  them  more  clearly. 
A  simple  explosive  is  composed  of  the  explosive  material  and  the  gas  produced  during 
the  reaction.  This  model  follows  the  basic  equations  of  conservation  for  inviscid  and 
adiabatic  flow.  The  conservation  equations  utilized  for  the  one-phase  1-D  model 
are:  Equation  10  conservation  of  mass,  Equation  11  conservation  of  momentum, 
and  Equation  12  conservation  of  energy  [5].  Mass,  momentum,  and  energy  is  not 
exchanged  between  the  gas/solid  mixture  and  the  environment.  Viscosity  and  thermal 
conductivity  can  be  ignored  because  it  occurs  over  a  longer  time  scale  than  that  of 
an  explosion,  which  occur  in  a  matter  of  microseconds.  Since  viscosity  and  thermal 
conductivity  are  ignored,  the  flow  can  be  described  by  the  Eulerian  Conservation 
Equations  10-12  below: 
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where  x  is  position,  t  is  time,  p  is  mixture  density,  u  is  velocity,  P  is  pressure,  and  E 
is  total  specihc  energy. 

In  order  apply  the  Euler  equations  to  this  system,  an  EOS  must  be  utilized  to 
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close  the  system.  Equation  13  below  is  a  caloric  EOS,  where  specihc  internal  energy 
is  a  function  of  both  density  and  pressure. 

e  =  e{p,P)  (13) 

where  e  is  the  mixture  specihc  internal  energy. 

3.2.2  Two-Phase  Model 

The  two-phase  model  is  similar  to  the  one-phase  model  in  that  both  are  for  simple 
explosive  and  use  the  conservation  equations,  bnt  now  both  phases  must  me  consid¬ 
ered.  For  the  two-phase  model,  a  simple  reactant/prodnct  mixtnre  composes  two 
phases,  a  solid  and  a  gas  phase.  The  conservation  equations  are  then  separated  for 
each  component.  The  volnme  fraction  (0)  of  each  component  is  inclnded  in  the  eqna- 
tions.  Thronghout  the  report  both  mass  fraction  (A)  and  0  are  ntilized  to  describe  the 
qnantity  of  alnminnm  within  the  high  explosive  material.  Both  can  be  calculated  by 
using  a  ratio  of  densities  of  the  materials  in  the  explosive  mixtnre.  BN  formnlated  a 
system  of  eqnations  for  basic  two-phase  flow  explosive  modeling,  shown  by  Equations 
14-20  [3]. 

The  terms  on  the  right  side  of  Eqnations  14-19  are  called  “source  terms”.  The 
mass  source  term  (C)  indicates  the  mass  exchange  between  the  solid  and  gas  phase, 
the  momentnm  sonrce  term  (Al)  indicates  the  momentnm  exchange  between  phases, 
and  the  energy  sonrce  term  (S)  indicates  the  energy  exchange  between  phases.  The 
conservation  of  mass  equation  for  the  solid  phase,  Eqnation  14,  has  a  positive  C  and 
the  conservation  of  mass  eqnation  for  the  gas  phase,  Eqnation  15,  has  a  negative  C. 
These  opposite  terms  between  the  solid  and  gas  phase  indicate  that  what  is  gained 
by  one  phase  must  be  lost  by  the  other  phase  or  vice  versa  to  maintain  conservation. 
The  relationship  between  phases  is  applied  to  the  six  conservation  eqnations;  Eqnation 


14-15,  Equation  16-17,  and  Equation  18-19. 
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where  x  is  position,  t  is  time,  pg  is  solid  volume  fraction,  pg  is  gas  volume  fraction, 
Ps  is  solid  density,  pg  is  gas  density,  Ug  is  solid  particle  velocity,  Ug  is  gas  particle 
velocity,  Pg  is  solid  pressure,  Pg  is  gas  pressure.  Eg  is  solid  total  specific  energy,  and 
Eg  is  gas  total  specific  energy.  The  subscripts  s  indicates  the  solid  components  and 
subscript  g  indicates  the  gas  components. 

Equation  20  relates  the  volume  fraction  exchange  in  the  reaction.  The  compaction 
source  term  (J^)  associates  the  inert  pore  collapse  of  the  granular  material  during  the 
reaction.  Equation  21  states  that  the  volume  fraction  of  solid  and  gas  components 
must  equal  to  one,  assuming  there  is  no  other  materials  or  massless  voids  in  the 
mixture. 
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Just  like  the  one-phase  model,  an  EOS  for  each  phase  is  needed  to  close  the  system 
of  Equations  14-20.  These  EOS  are  given  by  Equation  22  and  23. 

es  =  es{ps,Ps)  (22) 


~  ^giPgi  Pg)  (23) 

where  e*  is  the  solid  specihc  internal  energy  and  Cg  is  the  gas  specihc  internal  energy. 

In  order  to  develop  an  accurate  model  the  formulas  for  the  source  terms,  C,  M, 
and  T  must  be  determined.  The  formulations  for  the  source  terms  were  derived  from 
the  Second  Law  of  Thermodynamics  by  Bdzil  et  al.  [13] .  These  source  term  formulas 
are  shown  by  Equations  24-26.  In  these  source  term  equations  there  are  two  different 
parameters  of  interest.  The  relaxation  coefficients  and  the  partitioning  parameters. 

Relaxation  Coefficients 

The  relaxation  coefficients  “control  the  rate  at  which  the  solid  and  gas  pressures, 
velocities,  and  temperatures  equilibrate  [6].”  The  compaction  viscosity  ffisg)  is  uti¬ 
lized  to  determine  the  rate  at  which  the  mixture  mechanically  equilibrates,  the  drag 
coefficient  {Sgg)  is  utilized  to  govern  the  rate  of  kinematic  equilibration,  and  the 
heat  transfer  coefficient  {Hgg)  is  utilized  to  dictate  the  rate  of  thermal  equilibration 
[6].  In  a  piston  impact  test,  when  the  piston  hits  the  explosive  mixture  a  detonation 
wave  forms  and  travels  throughout  the  mixture.  At  any  position  in  the  mixture,  after 
the  detonation  wave  passes  the  mixture  will  relax  to  a  lower  state.  The  relaxation 
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coefficients  are  estimated  from  explosive  empirical  relaxation  data.  The  larger  these 
coefficients  are,  the  faster  the  mixture  will  equilibrate  after  the  wave  passes.  The 
smaller  these  coefficients  are,  the  longer  it  will  take  to  equilibrate. 

Partitioning  Parameters 

Partitioning  parameters  are  similar  to  the  relaxation  coefficients  in  the  way  that  they 
describe  the  interaction  between  phases  like  solid-gas.  The  partitioning  parameters 
determine  the  fraction  of  thermal  energy  exchange  that  occurs  during  different  pro¬ 
cesses  like  the  compaction,  chemical  reaction,  drag,  and  other  additional  dissipative 
processes  [6].  The  compaction  partitioning  function  {cgg)  determines  the  thermal 
energy  exchange  due  to  the  compaction  process,  the  chemical  reaction  partitioning 
function  {^sg)  determines  the  thermal  energy  exchange  due  the  chemical  reaction  pro¬ 
cess,  and  the  drag  partitioning  function  determines  the  energy  exchange  due  to 
the  drag  produced  by  the  production  of  combustion  gas  [6].  It  is  recognized  “that 
the  original  BN  model  implicitly  assigned  all  dissipation  energy  from  some  process 
to  either  the  gas  or  solid,  which  results  in  nonphysical  phase  temperature  increases. 
[6]”  Equations  24-26  are  the  formulations  for  the  source  terms  were  derived  from  the 
Second  Law  of  Thermodynamics  by  Bdzil  et  ah  [13]. 
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where  P  is  the  compaction  source  term,  A1  is  the  momentum  source  term,  C  is  the 
mass  source  term,  E  is  the  energy  source  term,  [igg  is  the  compaction  viscosity,  bgg  is 
the  drag  coefficient,  Plgg  is  the  heat  transfer  coefficient,  agg  is  the  drag  partitioning 
function,  Cgg  is  the  compaction  partitioning  function,  ^gg  is  the  chemical  reaction 
partitioning  function,  jigg  is  the  inter-granular  stress  from  particles,  Tg  is  the  solid 
temperature,  Tg  is  the  gas  temperature,  and  subscript  sg  indicates  the  relaxation 
process  between  the  solid  and  gas  phases. 


3.2.3  Multi-Phase  Model  with  N  Solid  Component 

The  multi-phase  model  is  derived  from  the  same  system  of  equations  as  the  two- 
phase  model  but  with  the  addition  of  more  equations  to  represent  the  other  added 
components.  The  model  provides  the  means  to  perform  simulations  with  mixtures 
with  multiple  solid  components.  A  multiple  solid  component  can  be  several  different 
combinations  of  solids.  One  combination  can  be  a  mixture  of  different  particle  sizes 
of  the  same  explosive.  Another  combination  can  be  a  mixture  of  different  explosive 
components.  The  combination  this  research  is  concerned  about  is  the  mixture  of 
explosive  and  metal  mixture.  The  conservation  equations  will  be  the  same  as  the 
two-phase  model  but  there  will  be  three  equations  for  the  gas:  conservation  of  mass, 
momentum,  and  energy.  There  will  also  be  three  conservation  equations  for  each  solid 
and  one  volume  fraction  equation  for  each  solid.  The  total  number  of  equations  need 
can  be  summarized  by  3-I-4N  equations,  where  N  is  the  number  of  solid  components 
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[6].  For  example,  this  research  will  focus  on  one  explosive  (PBX9501)  and  one  metal 
solid  {Al)  mixture  that  will  have  a  system  of  eleven  conservation  equations  because  of 
the  two  solids  (N=2).  Equations  27-33  were  derived  from  the  BN  formulated  system 
of  equations  for  this  multi-phase  model  with  N  number  of  solid  components  [3]. 
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where  (pi  is  i-th  solid  volume  fraction,  pi  is  i-th  solid  density,  Ui  is  i-th  solid  particle 
velocity.  Pi  is  i-th  solid  pressure,  and  Ei  is  i-th  solid  total  specihc  energy.  Other 
expressed  symbols  with  a  subscript  i  denote  parameters  of  the  solid  components. 

The  saturation  equation.  Equation  34,  shows  the  volume  fractions  for  all  the  i-th 
solids  and  the  volume  fraction  for  the  gas  component  equals  to  one. 
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The  “source  terms”,  C,M,  T,  for  all  the  components  are  summed  by  Equations  35- 
37.  These  equations  indicate  there  is  conservation  between  all  sources.  For  example, 
the  mass  sources,  C,  are  conserved  between  gas  and  all  solid  components. 


N 

Cg  +  Y.Ci  =  0  (35) 
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The  Tth  solid  mass  source  term  (Cj),  i-th  solid  moment  source  term  {Aii),  and 
i-th  solid  energy  source  term  (Si)  terms  can  be  separated  into  two  contributions; 
the  solid/gas  interaction  contributions  denoted  by  subscript  ig  and  the  solid/solid 
interaction  contribution  denoted  by  snbscript  im  [6].  Equations  38-40  below  show 
the  solid/gas  and  solid/solid  interaction  contributions  to  Ci,  Mi,  and 
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Ci  =  Sig  +  Sim  (40) 

m=l 

The  metal  ntilized  for  this  research  will  be  alnminum  and  will  be  treated  as  a 
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non-reactant  metal.  Since  non-reactant  Al  will  be  utilized  the  terms  that  include 
any  interaction  with  Al  will  not  be  utilized.  In  the  MPEXS  code  Al  is  denoted  by 
subscript  m  for  metal  and  the  explosive  is  denoted  by  s.  Applying  this  notation  to 
Equations  38-40  will  reduce  them  to: 

Cj  ^sg:  £i  Sgg 

3.2.4  Nozzling  Sources 

The  system  of  equations  formulated  by  BN  for  conservation  of  mass,  momentum, 
and  energy  of  explosive  reactions  can  be  written  in  a  general  form,  shown  by  Equation 
41  [6].  “The  mixture  entropy  inequality  gives  rise  to  terms  that  are  proportional  to 
the  volume  fraction  gradient  d(j)/dx  [6].”  These  non-conservative  terms  which  are 
multiplied  by  the  volume  fraction  gradients  are  known  as  nozzling  sources  (gj). 

i=l 

where  q  is  the  vector  of  unknown  variables  such  as  A,  tij,  Ej,  f  is  a  flux  vector  as 
a  function  of  q,  s  is  a  vector  of  algebraic  interphase  sources  as  a  function  of  q. 

The  left  side  of  Equation  41  are  the  Eulerian  conservation  equations  and  the  right 
side  are  the  source  term  equations.  “The  approximations  resulting  from  the  centered 
scheme  introduce  additional  numerical  error  into  the  solution.  This  is  due  to  the 
difficulties  associated  with  discretizing  the  nozzling  sources,  which  continues  to  be  an 
active  area  of  mathematical  research. [6]”  If  the  nozzling  terms  are  too  precise  then  it 
becomes  computationally  expensive,  if  they  they  are  not  precise  enough  then  it  will 
result  in  a  poor  approximation. 
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3.3  Multiphase  Explosive  Simulation  Code  Files 


3.3.1  Input  Files 

The  MPEXS  code  is  written  in  the  Fortran  programming  language.  The  code  is 
composed  of  several  routine  and  subroutine  hies.  The  purpose  of  the  MPEXS  code 
is  to  allow  the  user  to  input  the  explosive  simulation  parameters  in  as  few  hies  as 
possible.  The  latest  version  of  the  MPEXS  code  has  reduced  the  number  of  input 
hies  needed  to  run  a  simulation.  Most  of  the  inputs  are  done  in  the  maininput.txt 
hie  but  there  are  other  inputs  that  are  done  in  other  hies.  Brief  descriptions  are  given 
below. 


3.3. 1.1  maininput.txt 

The  maininput.txt  hie  is  is  the  principal  hie  where  most  of  the  simulation 
changes  will  be  made.  The  main  input  hie  is  divided  into  two  sections;  Simula¬ 
tion  Inputs  and  Initial  Conditions.  The  main  inputs  that  will  be  modihed  are  the 
number  of  solid  species  (Nsp),  number  of  discrete  computational  cells  (Njx),  left  and 
right  domain  boundary  locations  (LI)  and  (L2),  initial  reference  piston  speed  (vfO), 
hnal  simulation  time  (thnal),  average  component  grain  diameter  (dmean),  and  the 
types  of  equations  of  state  (eostypes).  The  initial  conditions  section  allows  the  user 
to  modify  the  volume  fraction  0  of  the  solids,  the  density  p  and  the  pressure  P  for 
both  the  gas  and  solids. 

The  MPEXS  code  can  run  both  neat  and  metalized  explosives.  The  term  “neat” 
is  utilized  to  refer  to  homogeneous  explosives  with  no  metal  particulates.  For  neat 
explosive  simulations  the  number  of  solid  species  Nsp  =  1.  By  doing  so,  the  code  will 
neglect  any  aluminum  properties  in  the  input  hies.  For  metalized  explosives  set  Nsp 
=  2. 

The  current  version  of  the  MPEXS  code  has  hve  diherent  equations  of  state;  Ideal 
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(0),  Virial  (1),  Mie-Griineisen  (2),  Wide- Ranging- Reactant  (4),  and  Wide- Ranging- 
Product  (5).  The  Ideal  EOS  are  applied  to  ideal  gases  in  many  flow  regimes.  Virial 
EOS  are  applied  to  non-ideal  gases  like  detonation  products  of  explosives.  Mie- 
Griineisen  EOS  applies  relationships  between  pressure  and  volume  of  solids  like  the 
metal  particulates. 

The  equations  of  state  can  be  selected  on  line  52  of  the  maininput.txt  hie.  The 
unreacted  explosive  PBX9501  will  use  the  Wide-Ranging-Reactant  equations  of  state, 
number  (4).  The  granular  aluminum  will  use  the  Mie-Griineisen  equations  of  state, 
number  (2).  The  gaseous  combustion  products  will  use  the  Wide- Ranging-Product 
equations  of  state,  number  (5).  When  running  cases  with  pure  PBX9501  explosive 
without  aluminum,  the  equations  of  state  are  called  in  line  52.  This  is  done  by 
setting  eostypes  =  4,  5,  0.  The  hrst  number  is  the  unreacted  explosive,  the  second 
number  is  the  gaseous  products,  and  the  third  number  is  a  place  holder.  The  third 
number  will  not  be  read  because  Nsp  is  set  to  1  for  pure  explosive  simulations.  When 
running  cases  with  PBX9501  and  aluminum,  line  ^  52  will  be  set  to  eostypes  =  4,  2, 
5.  The  hrst  number  is  the  unreacted  explosive,  the  second  number  is  the  aluminum, 
and  the  third  number  is  the  gaseous  products. 

3. 3. 1.2  timemarch.f90 

This  subroutine  steps  the  solution  forward  in  time  using  a  2nd  order  Strang  split¬ 
ting  method,  which  alternates  between  the  convective  and  local  source  portions  of 
the  solver  [6].  The  most  input  in  this  hie  is  the  frequency  data  is  stored.  The  default 
is  that  a  fort.xxx  hie  is  saved  every  10  time  steps.  This  input  is  found  in  line  ^  64 
that  reads  as  follows:  10)  ==  Q.OR.j  ==  Nts)THEN” .  An  analysis 

will  be  performed  to  determine  if  the  amount  of  data  collected  is  worth  the  increase 
in  accuracy. 
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3.3. 1.3  ppoutput.txt 


The  post  processing  output  file,  ppoutput.txt,  is  a  self  generated  file  needed  to 
calculate  the  detonation  wave  speed.  This  hie  is  stand-alone  and  can  be  modihed 
to  determine  the  wave  speed  at  different  times  intervals  by  changing  the  range  of 
fort  .XXX  hies.  If  the  MPEXS  code  stops  before  the  run  is  completed,  the  post 
processing  subroutine  will  not  generate  the  detonation  wave  speed  data.  In  that  case, 
the  inputs  can  be  manually  entered  and  be  able  to  calculate  detonation  wave  speeds. 
The  following  are  the  inputs;  number  of  discrete  grid  points,  number  of  condensed 
phases,  number  of  output  text  hies,  number  of  time  steps,  piston  velocity,  and  domain 
length. 


3.3. 1.4  postprocess. f90 

This  subroutine  performs  the  detonation  wave  speed  calculations  after  all  the  out¬ 
put  hies  have  been  created.  This  subroutine  is  can  also  be  utilized  as  a  stand-alone 
program  if  the  output  hies  have  already  been  created.  One  input  is  the  threshold  pres¬ 
sure.  This  threshold  pressure  is  utilized  to  diherentiate  between  a  passing  detonation 
wave  and  a  compaction  wave.  When  the  threshold  is  meet,  the  leading  wave  position 
at  each  recorded  time  will  be  stored  in  the  wavespeeddata.txt  hie.  The  threshold 
pressure  requires  some  trial  and  error  to  set  it  properly  for  detonation  waves.  If  too 
large  then  no  wave  will  be  detected  and  if  too  low  then  it  will  detect  any  pressure  rise 
even  if  it  is  not  a  detonation  wave.  Another  input  is  the  number  of  points  utilized  for 
the  linear  ht  that  determines  the  detonation  wave  speed.  The  more  points  utilized  the 
better  the  linear  ht.  If  too  many  points  are  use  the  result  will  be  inaccurate  because 
non  steady  state  points  will  be  utilized. 
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3.3.2  Execution  Files 


The  MPEXS  hydrocode  is  executed  by  using  the  Linux  command  terminal.  The 
binscript.sh  shell  script  is  utilized  to  execute  all  the  subroutines  and  input  hies.  As 
the  code  runs,  the  current  time  and  time  step  is  displayed  on  the  terminal  screen. 
The  data  displayed  on  the  terminal  screen  is  useful  to  determine  if  the  code  is  running 
properly.  The  data  is  saved  in  fort.xxx  hies  every  10  time  steps  which  is  the  default 
but  it  can  be  changed  under  the  timemarch.fOO  hie. 

3.3.3  Output  Files 

The  data  is  saved  to  output  hies  named  fort.xxx.  The“xxx”  is  a  number  sequence 
that  begins  at  100  and  increases  by  one  for  each  set  of  data.  Each  set  of  data  is 
collected  every  10  time  steps  and  it  depicts  what  is  happening  at  that  time  step. 
The  fort.  101  hie  will  contain  data  at  the  20th  time  step,  fort.  102  will  contain  data 
at  the  30th  time  step,  and  so  on.  The  fort.xxx  hies  contain  columns  of  data  for 
time,  position,  density,  pressure,  temperature,  volume  fraction,  and  particle  velocity. 
Table  1  shows  the  names  corresponding  to  each  column  of  data  for  a  neat  explosive 
simulation  which  includes  the  explosive  solid  and  gas  reaction  products.  Table  2 
shows  the  names  corresponding  to  each  column  of  data  for  a  metalized  simulation 
which  include  the  explosive  solid,  the  metal  solid,  and  the  gas  reaction  products. 

Table  1.  Fort.xxx  Columns  for  Two  Components  (Explosive  solid,  Gas  products) 
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where  t  is  time  in  microseconds,  x  is  position  in  millimeters,  p  is  density  in  kilo¬ 
grams  per  meters  cubed,  P  is  pressure  in  Pascals,  T  is  temperature  in  Kelvin,  0  is 
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volume  fraction,  and  u  is  particle  velocity  in  kilometers  per  second.  Subscript  “e” 
represents  explosive  solid  and  “g”  represents  gas  products. 

Table  2.  Fort.xxx  Columns  for  Three  Components  (Explosive  solid,  Metal  solid,  Gas 
products) 
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where  t  is  time  in  microseconds,  x  is  position  in  millimeters,  p  is  density  in  kilo¬ 
grams  per  meters  cubed,  P  is  pressure  in  Pascals,  T  is  temperature  in  Kelvin,  0  is 
volume  fraction,  and  u  is  particle  velocity  in  kilometers  per  second.  Subscript  “e” 
represents  explosive  solid,  “m”  represents  metal  solid,  and  “g”  represents  gas  prod¬ 
ucts. 

The  wavespeeddata.txt  hie  is  a  compilation  of  data  collected  from  all  the 
fort  .XXX  hies.  The  time  and  position  of  the  shock  wave  is  collected  from  each 
fort  .XXX  hie  by  setting  a  threshold  pressure  to  detect  when  the  shock  wave  has  passed. 
The  time  and  position  are  saved  into  a  separate  hie  named  wavespeeddata.txt.  Ta¬ 
ble  3  below  provides  an  example  of  the  data  contained  in  the  wavespeeddata.txt 
hie.  The  data  is  stored  in  four  columns.  The  hrst  column  it  the  time  when  the  shock 
wave  passes,  the  second  column  is  the  position  of  the  shock  wave  in  the  explosive  par¬ 
ticles,  the  third  column  is  the  position  for  the  metal  particles,  and  the  forth  column 
is  the  position  for  the  gas  particles. 


Table  3.  Wavespeeddata.txt  Example 
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In  some  cases  a  detonation  will  not  occur  and  only  the  wave  due  to  the  impact 
of  the  piston  will  go  through.  This  wave  is  called  the  compaction  wave.  Detonation 
wave  velocities  are  much  larger  than  compaction  wave  velocities  and  will  be  noticeable 
when  a  detonation  does  not  occur.  The  goal  is  to  choose  a  velocity  threshold  large 
enough  that  will  detect  the  detonation  wave  but  ignore  the  compaction  wave.  The 
wavespeeddata.txt  data  provides  the  time  and  position  of  the  wave.  When  plotted 
the  detonation  wave  velocity  for  the  explosion  can  calculated  with  the  slope  of  the 
line  dx/dt.  Since  the  piston  face  is  utilized  as  the  reference,  the  piston  velocity  must 
be  added  to  the  wave  velocity.  The  piston  and  wave  velocities  will  give  the  detonation 
velocity. 

3.4  Simulation  Run  Setup 

In  order  to  study  the  effects  of  aluminum  particulates  on  the  detonation  veloc¬ 
ity  of  explosives,  computational  simulations  must  be  performed  with  various  particle 
size  and  mass  fractions.  This  research  will  focus  on  the  high  explosive  compound 
PBX9501.  The  hrst  step  is  to  test  the  hdelity  of  the  MPEXS  code  on  a  simple  ex¬ 
plosive  like  PBX9501  with  no  aluminum  particulates,  this  will  be  referred  as  neat 
PBX9501.  If  the  simulation  results  for  neat  PBX9501  have  a  good  correlation  with 
empirical  results,  then  the  second  step  is  to  run  simulations  with  aluminum  particu¬ 
lates.  Comparing  simulation  data  to  empirical  data  does  not  provide  any  beneht  to 
this  research  due  to  the  complexity  of  plastically  bonded  explosives  like  PBX9501. 
The  detonation  properties  are  different  for  every  batch  made  due  to  the  purity  and 
curing  process  of  the  binding  material  can  vary  from  batch  to  batch.  Instead  of  com¬ 
paring  simulation  data  to  empirical  data,  this  research  will  look  at  the  stability  of  the 
code  by  running  several  tests  to  determine  consistency  and  limits  of  the  code.  This 
research  will  look  at  trends  in  several  detonation  properties  for  metalized  PBX9501. 
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AFRL  Munitions  Directorate  currently  uses  spherical  aluminum  particulates  man¬ 
ufactured  by  Toyal  America  Inc  [12]  in  explosive  mixtures  being  tested.  Toyal  manu¬ 
factures  a  range  of  aluminum  powers  in  different  diameters  for  aerospace  and  military 
applications.  Toyal  manufactures  six  different  sized  aluminum  powders  with  a  mean 
diameter  ranging  from  9  to  34  microns  [12].  The  six  different  diameters  will  be 
utilized  as  a  varying  parameter  for  the  computational  simulations  of  this  research. 

Another  parameter  that  will  be  varied  is  the  mass  fraction  (A)  of  aluminum  to 
the  solid  explosive.  Increasing  the  mass  fraction  of  the  aluminum  will  decrease  the 
mass  fraction  of  the  solid  explosive.  A  decrease  in  the  amount  of  explosive  available 
for  the  reaction  will  affect  detonation  properties  such  as  detonation  velocity,  run-to- 
detonation  time,  and  run-to-detonation  location.  Three  different  mass  fractions  will 
be  studied  for  this  research;  10,  15,  and  20  percent. 

Figure  5  below  shows  the  two  varying  parameters  that  will  be  studied,  aluminum 
mass  fraction  and  aluminum  mean  diameter.  The  first  set  of  simulations  that  will 
be  performed  will  be  with  neat  PBX9501.  The  data  from  these  simulations  will 
utilized  verify  that  the  equations  of  state  in  the  MPEXS  hydrocode  are  correct.  The 
verihcation  of  the  MPEXS  code  can  be  done  by  comparing  run-to-detonation  results  to 
other  hydrocodes.  The  next  set  of  simulations  will  consider  different  cases  of  PBX9501 
with  aluminum.  The  six  aluminum  diameters  will  be  tested  for  each  aluminum  mass 
fraction.  For  example,  the  aluminum  mass  fraction  of  10  percent  will  be  tested  with 
9,  13,  19,  23,  29,  and  34  micron  diameter.  The  same  simulations  will  be  performed 
for  15  and  20  percent. 
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Neat  PBX9501 
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Figure  5.  Varying  Parameters  for  Simulation  Runs 


3.5  Neat  Explosive  Simulations 

Two  things  must  be  considered  before  the  study  of  the  effects  of  metalized  ex¬ 
plosives;  the  validity  of  the  MPEXS  code  and  a  baseline  of  initial  parameters  to  be 
utilized.  In  this  section,  different  computational  cell  quantities  and  different  piston 
speeds  will  studied  with  the  purpose  to  determine  a  baseline  of  initial  parameters  to 
be  utilized  for  metalized  explosives. 

For  the  neat  explosive  simulations,  data  will  be  collected  at  different  piston  speeds 
{up)  and  at  different  number  of  discrete  computational  cells  {Nx)  for  both  PBX9501 
and  HMX  explosives.  Simulations  with  small  number  of  computational  cells  will 
run  in  a  shorter  time  but  the  resolution  will  be  lower.  The  lower  resolution  might 
miss  important  data  points  specially  at  the  location  of  the  shock  wave.  Simulations 
with  very  large  number  of  computational  cells  will  provide  a  higher  resolution  with 
more  data  points  around  the  shock  but  will  take  longer  to  run.  Besides  run  time 
and  resolution,  the  amount  of  data  output  per  simulation  is  an  another  important 
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consideration.  The  more  computational  cells  the  more  output  hies  the  code  will 
produce.  Simulation  run  time  is  not  as  important,  but  data  output  is  important. 
Therefore,  a  study  will  be  performed  to  determine  if  the  higher  resolution  simulations 
are  worth  the  data  each  simulation  will  output  based  on  the  percent  difference  for 
each  number  of  computational  cells.  Detonation  properties  can  vary  depending  on 
the  impact  speed  of  the  piston.  A  study  will  be  performed  to  determine  what  piston 
speed  will  provide  the  best  results  and  what  piston  speeds  will  cause  the  code  to  fail. 

For  the  two  initial  studies,  computational  cells  and  piston  speed,  a  low,  medium, 
and  higher  limit  will  be  tested  to  gauge  the  detonation  properties  of  neat  explosives. 
For  computational  cell,  runs  with  800  and  2000  cells  will  be  the  lower  and  higher 
limits  to  be  tested.  For  piston  speed,  runs  with  0.5,  1.0,  2.0  km/s  will  be  tested. 
Both  the  computational  cell  quantity  and  piston  speed  studies  are  described  in  more 
detail  in  the  following  sections. 

3.5.1  PBX9501  at  Different  Computational  Cell  Quantities 

Two  different  number  of  discrete  computational  cells  will  be  considered  for  each 
explosive,  800  and  2000  cells.  The  purpose  is  to  determine  the  amount  of  cells  needed 
to  obtain  good  results  and  if  the  run  time  and  data  storage  of  the  simulation  is  worth 
the  resolution.  In  this  case,  the  explosive  PBX9501  was  simulated  at  three  different 
piston  speeds  and  the  two  computational  cell  quantities  with  the  focus  of  comparing 
results  for  800  and  2000  cells.  Larger  numbers  of  computational  cell  sizes  such  as  2400, 
2500,  3000  cells  were  tested  but  the  code  did  not  run  with  the  current  parameters.  It 
is  believed  that  computational  limits  on  the  computer  were  reached  and  the  computer 
aborts  the  code.  This  issue  will  be  discussed  in  Results  chapter  after  the  code  is  run 
on  the  Air  Force  Institute  of  Technology  (AFIT)  computer  cluster  and  the  results  are 
analyzed.  There  are  hve  detonation  properties  that  are  considered  through  out  this 
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research;  density,  pressure,  temperature,  particle  velocity,  and  volume  fraction  of  the 
solids  and  gas. 

When  the  simulation  was  run  at  a  piston  speed  of  Up  =  0.5  km/s  the  results 
did  not  seemed  to  vary  from  800  to  2000  computational  cells  for  density,  pressure, 
temperature,  particle  velocity,  and  volume  fraction.  Simulations  with  800  cells  took 
around  13  minutes  to  run  and  stored  about  746  MB.  Simulations  with  2000  cells  took 
around  78  minutes  to  run  and  stored  about  4.6  gibibyte  {GiB).  As  a  reference,  1.0 
gibibyte  equals  1.07374  gigabytes.  These  run  times  and  storage  sizes  are  for  neat 
explosive,  where  no  Al  is  present.  Figure  6a  shows  pressure  versus  position  for  800 
computational  cells  and  Figure  6b  shows  the  same  plot  for  2000  computational  cells. 
As  seen  in  the  hgures  the  plotted  data  looks  identical.  When  the  data  is  compared, 
it  can  be  seen  that  there  is  a  1.5%  difference  between  800  cells  and  2000  cells.  The 
number  of  computational  cells  utilized  for  0.5  km/s  makes  a  very  small  difference 
in  the  results.  A  generalized  conclusion  about  the  effects  of  computational  cells  can 
not  be  made  at  this  point  based  on  the  results  because  the  detonation  did  not  reach 
steady  state.  There  are  two  options,  either  run  the  simulation  for  a  longer  time  or 
increase  the  piston  impact  speed.  To  keep  the  position  domain  consistent,  faster 
piston  speeds  were  considered  and  are  discussed  in  the  next  paragraph. 
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(b)  PBX9501  WpO.5  2000cells 


Figure  6.  PBX9501  P  vs  X  Piston  Speed  =  0.5  km/s 


Since  steady  state  was  not  reached  within  the  position  domain  for  piston  speed 
of  0.5  km/s,  the  piston  speed  was  increased.  Simulations  with  piston  speeds  of  1.0 
and  2.0  km/s  were  performed  for  PBX9501.  The  focus  is  on  the  grid  resolution 
therefore  the  same  computational  cells,  800  and  2000,  were  utilized  for  the  faster 
piston  speeds.  800  and  2000  computational  cells  were  utilized  for  the  simulations 
with  faster  piston  speeds.  Since  the  results  for  1.0  km/s  and  2.0  km/s  are  similar 
when  varying  computational  cell  sizes,  only  the  2.0  km/s  case  will  be  presented  and 
discussed  here.  The  results  for  density,  pressure,  particle  velocity,  and  volume  fraction 
were  almost  identical  when  800  and  2000  computational  cells  are  compared.  For  1.0 
and  2.0  km/s  the  detonation  reached  a  steady  state.  Figure  7a  and  Figure  7b  show 
pressure  vs  position  for  both  cases.  Comparing  the  results  for  the  two  cell  sizes  it  is 
calculated  there  is  a  1.0%  difference  in  pressure  results,  0.8%  difference  for  particle 
velocity,  1.5%  difference  for  density,  and  0.0%  difference  for  detonation  velocity.  The 
difference  is  small  enough  that  it  can  be  concluded  that  computational  cells  sizes  do 
not  affect  these  detonation  parameters. 
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Pressure  vs  Position 


Pressure  vs  Position 


(a)  9501  Up2.0  SOOcells  (b)  9501  Up2.0  2000cells 

Figure  7.  9501  P  vs  X  Piston  Speed  =  2.0  km/s 

As  seen  in  Figure  8a  and  8b,  the  temperature  of  the  explosive  particles  has  a  higher 
spike  when  2000  computational  cells  are  utilized  versus  the  800  cells.  The  thin  region 
behind  the  detonation  wave  is  typically  very  small  compared  to  the  grid  sizes  utilized. 
The  more  computational  cells  utilized,  the  finer  the  grid  size  gets  and  more  data  can 
be  captured  in  this  region.  The  case  with  2000  cells  provides  more  detail  in  this  region 
and  shows  that  the  temperature  is  actually  higher.  When  less  computational  cells  are 
utilized,  some  of  the  higher  temperature  data  points  are  skipped,  therefore  showing 
a  lower  temperature.  Another  observation  for  temperature  vs  position  is  that  the 
transition  to  equilibrium  after  the  shock  passes  is  smoother  with  2000  computational 
cells.  The  smoother  transition  is  also  due  to  the  increase  of  data  points  per  unit 
distance. 
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Temperature  vs  Position 


Temperature  vs  Position 


(a)  9501  Up2.0  SOOcells 


(b)  9501  Up2.0  2000cells 


Figure  8.  9501  T  vs  X  Piston  Speed  =  2.0  km/s 


The  results  showed  there  was  repeated  trends  and  consistent  results  between  the 
two  different  computational  cell  sizes.  There  were  no  abnormal  results  between  the 
two  cell  sizes.  These  results  help  validate  that  the  code  is  providing  consistent  and 
predictable  results. 

3.5.2  HMX  at  Different  Computational  Cell  Quantities 

Neat  explosive  HMX  was  also  studied  at  different  piston  speeds  and  two  different 
computational  cell  quantities,  800  and  2000  cells.  To  compare  both  the  computational 
cell  quantities  a  piston  speed  of  0.5  km/s  was  considered  hrst.  The  simulation  of 
neat  HMX  with  a  piston  speed  of  0.5  km/s  showed  a  2.2%  difference  between  800 
and  2000  cells  in  particle  velocities,  0.7%  difference  in  temperatures,  1.7%  difference 
in  pressures,  0.3%  difference  in  densities,  and  2.3%  difference  in  detonation  velocities. 
Figure  9  shows  a  comparison  between  800  and  2000  cells  for  temperature  vs  position 
of  HMX.  Both  the  HMX  and  PBX9501  simulations  at  0.5  km/s  showed  that  there 
is  a  very  small  difference  between  800  and  2000  computational  cells.  Based  on  the 
results  for  a  piston  speed  of  0.5  km/s,  the  detonation  did  not  reach  steady  state. 
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No  conclusions  can  be  made  from  these  results  about  the  different  cell  quantities. 
To  make  a  good  comparison  of  computational  cell  sizes,  higher  piston  speeds  were 
considered. 


Temperature  vs  Position 


Temperature  vs  Position 


(a)  HMX  UpO.5  SOOcells 


(b)  HMX  UpO.5  2000cells 


Figure  9.  HMX  T  vs  X  Piston  Speed  =  0.5  km/s 


The  results  of  the  simulations  for  neat  HMX  with  piston  speed  of  2.0  km/s  had  less 
than  2%  difference  from  800  to  2000  computational  cells  for  density,  pressure,  particle 
velocity,  and  volume  fraction.  The  only  parameter  that  displayed  a  difference  between 
800  to  2000  cells  was  temperature  as  seen  in  Figure  10.  For  the  2000  computational 
cell  simulation,  the  temperature  reached  a  higher  temperature  and  had  a  smoother 
transition  behind  the  shock.  Like  the  PBX9501  simulation,  the  smoother  transition 
and  higher  temperature  is  due  to  the  larger  number  of  data  points  per  unit  length 
when  2000  cells  are  utilized. 
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3.5 


Temperature  vs  Position 


Temperature  vs  Position 


(a)  HMX  Up2.0  SOOcells 


(b)  HMX  Up2.0  2000cells 


Figure  10.  HMX  T  vs  X  Piston  Speed  =  2.0  km/s 


Based  on  this  study  for  neat  explosive  simulations,  it  is  shown  that  2000  cells  will 
increase  the  resolution  by  adding  more  data  points,  but  it  will  require  more  data  stor¬ 
age.  800  cells  provided  less  precise  results  but  it  required  less  storage.  The  difference 
in  the  results  between  800  and  2000  cells  is  within  1%  to  2%.  Running  the  simulations 
with  800  cells  is  a  better  option  for  faster  and  quick  checks  of  detonation  properties 
because  it  will  require  less  data  storage.  The  larger  data  storage  of  2000  cells  might 
be  an  issue  with  the  metalized  explosive  simulations  but  with  neat  explosives  is  not 
a  issue.  Therefore,  2000  cells  will  be  utilized  for  studying  the  different  speeds  of  neat 
explosives  in  the  next  section. 

3.5.3  PBX9501  at  3  Different  Piston  Speeds 

From  the  previous  study  it  was  determined  that  a  good  starting  resolution  was 
2000  computational  cells.  The  following  simulations  will  look  at  the  effect  of  different 
piston  speeds  to  determine  the  best  piston  speed  for  the  metalized  explosive  simula¬ 
tions.  It  is  important  to  get  a  good  understanding  of  the  neat  explosive  detonation 
properties  before  the  addition  of  metal  particulates. 
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Three  different  piston  speeds  were  chosen  for  the  following  simulations,  all  with 
2000  computational  cells.  Piston  speeds  of  0.5  km/s,  1.0  km/ s,  and  2.0  km/s.  All  hve 
detonation  properties  were  compared  at  each  of  the  piston  speeds.  The  detonation 
properties  are:  density,  pressure,  temperature,  particle  velocity,  and  volume  fraction. 

Figure  11  plots  the  volume  fraction  vs  position  for  the  three  different  piston  speeds. 
Figure  11a  shows  that  volume  fraction  does  not  transition  from  explosive  solid  to  gas 
for  the  Up  =  0.5  case.  The  explosive  volume  fraction  should  exponentially  decay  to 
zero  and  the  gas  volume  fraction  should  exponentially  increase  to  the  total  solids 
volume  fraction  as  shown  in  Figures  11b  and  11c.  As  discussed  earlier,  for  piston 
speed  of  0.5  km/s  the  detonation  wave  is  developing  and  does  not  reach  steady  state 
for  the  time  and  position  domain  considered.  Figure  11c  is  a  good  example  of  a 
detonation  wave  that  has  reached  its  steady  state.  Figure  11b  is  almost  at  steady 
state  but  not  quite  there.  Steady  state  is  reached  when  the  explosive  solid  decays  to 
zero  very  sharply  and  does  not  change. 
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Volume  Fraction  vs  Position 


(a)  PBX9501  UpO.5  2000grid 


Volume  Fraction  vs  Position  Volume  Fraction  vs  Position 


(b)  PBX9501  Mpl.O  2000grid  (c)  PBX9501  Up2.0  2000grid 

Figure  11.  PBX9501  Volume  Fractions 

The  rise  in  density  for  the  Up  =  0.5  km/s  case  is  due  to  the  compression  wave 
caused  by  the  impact  of  the  piston.  The  small  rise  in  density  indicates  the  piston 
speed  was  not  fast  enough  to  reach  steady  state  within  the  time  and  position  domain. 

Since  steady  state  was  not  reached  at  Up  =  0.5  km/s,  the  piston  speed  was  in¬ 
creased  to  1.0  km/s  and  2.0  km/s.  The  results  for  Up  =  1.0  km/s  and  2.0  km/s  are 
shown  in  Figure  12b  and  12c.  For  these  faster  speeds  it  can  be  seen  that  a  detonation 
occurred  because  there  is  a  larger  increase  in  density  as  the  shock  wave  passed.  The 
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peak  density  for  the  gas  is  2437  kg/m^  at  Up  =  1.0  km/ s  and  2484  kg/m^  at  Up  = 
2.0  km/s.  The  biggest  difference  between  the  two  speed  is  the  rate  the  rarefaction 
wave  reaches  eqnilibrinm.  The  Up  =  2.0  km/s  simnlation  reaches  eqnilibrinm  at  a 
faster  rate  and  the  Up  =  1.0  km/s  simnlation  at  a  slower  rate.  Fignre  12b  shows  a 
smoother  more  gradnal  decrease  to  the  eqnilibrinm  state  and  Fignre  12c  a  sharper 
drop  to  eqnilibrinm  state. 
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(a)  Density  vs  X  Up  =  0.5 


Density  vs  Position 


Density  vs  Position 


(b)  Density  vs  X  Up  =  1.0 


(c)  Density  vs  X  Up  =  2.0 


Figure  12.  PBX9501  Density  vs  X  3  Speeds 


Pressnre  versns  position  is  plotted  for  the  three  piston  speeds  and  the  resnlts  are 
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shown  in  Figure  13.  The  same  patterns  as  density  are  seen  for  pressure,  were  steady 
state  is  not  reached  at  Up  =  0.5  km/s  and  a  steady  state  detonation  occurs  for  the 
faster  piston  speeds.  The  peak  pressure  for  the  gas  is  38.86  GPa  at  Up  =  1.0  km/s  and 
42.64  GPa  at  Up  =  2.0  km/s.  Piston  speed  directly  affects  the  detonation  properties 
like  pressure.  When  the  piston  speed  was  increased,  the  pressures  also  increased. 
In  Figure  13b,  the  rarefaction  wave  for  Up  =  1.0  km/s  has  a  more  gradual  drop  in 
pressure.  Figure  13c  shows  the  rarefaction  wave  for  Up  =  2.0  km/s  has  a  sharper 
drop  in  pressure  and  then  equilibrates. 


Pressure  vs  Position 


(a)  Pressure  vs  X  Up  =  0.5 


Pressure  vs  Position  Pressure  vs  Position 


(b)  Pressure  vs  X  Up  =  1.0 


(c)  Pressure  vs  X  Up  =  2.0 


Figure  13.  PBX9501  Pressure  vs  X  3  Speeds 
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The  relationship  between  the  three  speeds  is  the  same  the  remaining  properties 
like  temperatnre,  particle  velocity,  and  volnme  fraction.  The  correlation  between 
piston  speed  and  detonation  properties  is  that  as  piston  impact  speed  increases,  the 
detonation  properties  increase.  Higher  piston  speeds  resulted  in  higher  pressure, 
temperature,  particle  velocity,  and  detonation  velocity.  For  brevity  of  the  report 
some  of  the  key  detonation  property  values  for  PBX9501  at  the  three  piston  speed 
are  shown  in  Table  4. 


Table  4.  Neat  PBX9501  Detonation  Properties  for  Different  Piston  Speeds 


Piston  Speed, 

Up  {km/s) 

Pressure, 
Pg  {GPa) 

Temperature, 
Tg  {K) 

Part.  Velocity, 

Ug  {km/s) 

Det.  Velocity, 

D  {km/s) 

0.5 

2.07 

969 

0.576 

2.252 

1.0 

38.86 

29170 

3.357 

7.842 

2.0 

42.64 

33280 

3.594 

7.842 

3.5.4  HMX  at  3  Different  Piston  Speeds 

The  same  analysis  was  performed  for  the  neat  HMX  at  three  different  piston 
speeds;  0.5  km/s,  1.0  km/s,  and  2.0  km/s.  When  the  piston  speed  was  increased, 
neat  HMX  had  the  detonation  properties  trends  as  neat  PBX9501.  As  the  piston 
speed  increased  the  pressure,  temperature,  particle  velocity,  and  detonation  velocity 
also  increased.  A  summary  of  the  results  for  neat  HMX  is  on  Table  5. 


Table  5.  Neat  HMX  Detonation  Properties  for  Different  Piston  Speeds 


Piston  Speed, 

Up  {km/s) 

Pressure, 
Pg  {GPa) 

Temperature, 
Tg  {K) 

Part.  Velocity, 
Ug  {km/s) 

Det.  Velocity, 

D  {km/s) 

0.5 

2.35 

1115 

0.585 

2.629 

1.0 

33.38 

22740 

2.716 

7.799 

2.0 

36.76 

27760 

2.869 

7.968 
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3.5.5  Neat  PBX9501  and  Neat  HMX  Detonation  Properties  Compar¬ 


ison 

The  previous  two  sections  looked  at  both  neat  PBX9501  and  neat  HMX  individu¬ 
ally  to  determine  how  computational  cells  and  piston  speed  affected  detonation  prop¬ 
erties  for  each  explosive.  This  section  will  focus  on  comparing  detonation  properties 
between  PBX9501  and  HMX  at  a  piston  speed  of  2.0  km/s  and  2000  computational 
cells.  If  the  code  is  running  correctly  then  the  detonation  properties  for  PBX9501 
and  HMX  should  have  similar  trends. 

Figure  14  below  plots  the  the  density  vs  the  position  for  both  PBX9501  and  HMX 
at  piston  speed  of  2.0  km/s  and  2000  computational  cells.  The  density  of  the  explosive 
material  as  the  detonation  wave  passes  is  higher  for  PBX9501.  The  stable  density 
values  of  the  rarefaction  wave  for  PBX9501  are  higher  as  well. 


Density  vs  Position 


Density  vs  Position 


(a)  HMX  Up2.0  2000grid 


(b)  9501  Up2.0  2000grid 


Figure  14.  HMX  and  9501  p  vs  X  Piston  Speed  =  2.0  km/s 


Figure  15  plots  the  pressure  vs  position  of  both  HMX  and  PBX9501.  The  pressure 
for  PBX9501  is  higher  than  the  pressure  for  HMX.  The  peak  pressure  for  PBX9501 
is  higher  than  that  of  HMX  at  piston  speed  of  2.0  km/s,  2000  cells. 
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Pressure  vs  Position  Pressure  vs  Position 


(a)  HMX  Up2.0  2000grid  (b)  9501  Up2.0  2000grid 

Figure  15.  HMX  and  9501  P  vs  X  Piston  Speed  =  2.0  km/s 

Figure  16  plots  the  temperature  vs  position  of  both  HMX  and  PBX9501.  The 
temperature  for  PBX9501  is  higher  than  the  temperature  for  HMX.  The  peak  tem¬ 
perature  for  PBX9501  is  higher  than  that  of  HMX  at  piston  speed  of  2.0  km/s,  2000 

cells. 


(a)  HMX  Up2.Q  2000grid 


(b)  9501  Up2.0  2000grid 


Figure  16.  HMX  and  9501  T  vs  X  Piston  Speed  =  2.0  km/s 


Figure  17  plots  the  velocity  vs  position  of  both  HMX  and  PBX9501.  The  velocity 
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for  PBX9501  is  higher  than  for  HMX.  The  peak  particle  velocity  for  PBX9501  is 
higher  than  that  of  HMX  at  piston  speed  of  2.0  km/ s,  2000  cells. 


Particle  Velocity  vs  Position  Particle  Velocity  vs  Position 


(a)  HMX  Up2.0  2000grid  (b)  9501  Up2.Q  2000grid 

Figure  17.  HMX  and  9501  u  vs  X  Piston  Speed  =  2.0  km/s 

Figure  18  plots  the  volume  fraction  vs  position  of  both  HMX  and  PBX9501.  The 
volume  fraction  for  both  HMX  and  PBX9501  are  very  similiar  if  not  the  same.  The 
detonation  for  both  explosives  are  at  steady-state  and  the  volume  fractions  should  be 
at  zero  for  the  explosive  solid  and  at  one  for  the  gas  products  after  the  detonation 
wave  has  passed. 
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Volume  Fraction  vs  Position  Volume  Fraction  vs  Position 


(a)  HMX  Up2.0  2000grid  (b)  9501  Up2.0  2000grid 

Figure  18.  HMX  and  9501  (/)  vs  X  Piston  Speed  =  2.0  km/s 

Neat  PBX9501  explosive  is  polymer  bonded  explosive  that  is  composed  of  95% 
HMX  and  bonnded  together  by  2.5%  estane  and  other  materials.  The  addition  of 
estane-based  binder  material  in  PBX9501  adds  more  fnel  to  the  chemical  reaction 
and  can  be  clearly  seen  in  the  stoichiometry  eqnations.  The  stoichiometry  calcnla- 
tions  show  that  more  fnel  is  available  to  the  reaction  dne  to  the  binder  materials. 
Detonation  properties  snch  as  pressnre,  temperature,  density,  and  particle  velocity 
are  all  affected  by  the  addition  of  the  binder  material.  This  effect  was  seen  in  Figure 
14  to  Figure  18.  The  binder  material  retards  the  burn  rate  of  the  explosive  and  this 
retardation  makes  the  reaction  last  longer  in  which  increases  the  detonation  properties 
like  density,  pressure,  temperature,  and  particle  velocity. 

3.6  Metalized  Explosive  Simulations 

3.6.1  PBX9501/A1  at  Different  Computational  Cell  Quantities 

A  few  cases  with  metalized  PBX9501  explosive  were  done  to  determine  the  accu¬ 
racy  of  the  higher  resolution  simulations  and  how  much  data  storage  is  needed  for 
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each  metalized  explosive  simulation.  Simulations  with  PBX9501  and  aluminum  were 
performed  at  a  piston  speed  of  1.5  km/s,  10%  Al  mass  fraction,  and  9/im  and  Sifim 
Al  particle  diameter.  Also,  simulations  at  Up  =  1.5  km/s,  20%  Al  mass  fraction,  and 
9pm  Al  particle  diameter.  These  three  simulations  will  cover  a  good  range  of  the 
cases  that  will  be  run  for  the  results  chapter. 

Figure  19  below  plots  pressure  vs  position  for  10%  Al  mass  fraction  and  9pm  to 
compare  800  and  2000  cell  simulations.  Figure  20  below  plots  temperature  vs  position 
for  10%  Al  mass  fraction  at  9pm  to  compare  800  and  2000  cell  simulations.  As  seen 
in  Figure  20,  the  gas  temperatures  rapidly  increase  to  values  outside  the  chart  range. 
This  large  spike  in  gas  temperature  are  consequence  of  the  computational  model 
utilized  for  the  simulations.  The  large  spike  were  also  seen  by  other  researchers  like 
Schwendeman,  Wahle,  and  Kapila  as  described  in  Dr.  Crochet’s  dissertation  [5]. 


Pressure  vs  Position 


Pressure  vs  Position 


(a)  PBX9501  10%  9fim  Upl.5  SOOcells  (b)  PBX9501  10%  9^m  Upl.5  2000cells 


Figure  19.  PBX9501  10%  9/im  P  vs  X  Piston  Speed  =  1.5  km/s 
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Temperature  vs  Position 


Temperature  vs  Position 


(a)  PBX9501  10%  9fim  Upl.5  SOOcells  (b)  PBX9501  10%  9^m  Upl.5  2000cells 

Figure  20.  PBX9501  10%  9^m  T  vs  X  Piston  Speed  =  1.5  km/s 

A  summary  of  the  simulations  performed  to  study  the  precision  of  higher  com¬ 
putational  cells  and  the  amount  of  storage  needed  per  simulation  are  listed  below  in 
Table  6  and  Table  7.  To  describe  each  simulation  in  a  shorter  format,  the  following 
naming  convention  is  utilized,  “PBX9501_10_9”.  The  hrst  part  is  the  name  of  the 
explosive,  the  hrst  number  after  the  explosive  name  is  the  mass  fraction  of  the  alu¬ 
minum  in  percent,  and  the  last  number  is  the  diameter  of  the  aluminum  in  microns. 
Other  added  descriptions  are  added  to  the  name.  The  added  description  of  “upl.5” 
is  the  piston  speed  and  “SOOcells”  is  the  number  of  computational  cells  utilized. 

Table  6  shows  the  percent  difference  between  each  simulation  at  800  and  2000 
computational  cells.  For  the  hrst  simulation  listed,  PBX9501_10_9_upl.5,  there  is 
a  1.5%  average  diherence  between  800  and  2000  cells  for  particle  velocities,  23.4% 
diherence  for  temperature,  2.1%  for  pressure,  0.6%  for  density,  and  0.5%  for  detona¬ 
tion  velocity.  The  other  two  cases  are  listed  in  the  table.  The  largest  diherence  is  the 
temperature  but  as  described  in  the  previous  paragraph,  the  computational  model 
does  not  accurately  predict  the  actual  temperatures  near  the  shock  wave.  With  this 
in  mind,  the  focus  will  be  on  the  other  detonation  properties.  The  other  properties 
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for  PBX9501_10_9_upl.5,  have  percent  differences  below  2.1%.  Besides  the  per¬ 
cent  difference  for  temperature,  PBX9501_20_9_upl.5  has  a  large  percent  difference 
in  detonation  velocity.  The  reason  for  the  large  percent  difference  compared  to  the 
other  cases  is  because  the  PBX9501_10_9_upl.5  800  cell  simulation  ended  earlier 
than  the  2000  cell  simulation.  The  next  paragraph  will  look  at  how  these  differences 
compare  to  the  amount  of  storage  needed  per  simulation. 


Table  6.  Difference  Between  800  and  2000  Computational  Cells 


Part.  Vel. 
%  Diff. 

Temp. 
%  Diff. 

Press. 

%  Diff. 

Density 
%  Diff. 

Det.  Vel. 
%  Diff. 

PBX9501_10_9_upl.5 

1.5% 

23.4% 

2.1% 

0.6% 

0.5% 

PBX9501_10_34_upl.5 

1.4% 

8.6% 

1.9% 

0.5% 

0.5% 

PBX9501_20_9_upl.5 

2.3% 

16.0% 

1.7% 

0.4% 

5.1% 

Table  7  below  shows  a  comparison  between  the  amount  of  storage  needed  for 
the  data  of  the  simulation  for  800  and  2000  computational  cells.  The  hrst  case 
PBX9501_10_9_upl.5  at  800  cells  outputted  2.9  gibibytes  and  2000  cells  outputted 
23.7  gibibytes.  The  2000  cell  simulation  requires  8.2  times  more  data  storage  capacity 
than  the  800  cell  simulation.  The  2000  cell  simulation  for  PBX9501_10_34_upl.5 
requires  6.5  times  more  data  storage  than  800  cells.  The  2000  cell  simulation  for 
PBX9501  20  9  upl.5  requires  7.7  times  more  data  storage  than  800  cells. 


Table  7.  Difference  Between  800  and  2000  Computational  Cells 


Data 

Output(GiB) 

2000  vs  800 
(Times  Larger) 

PBX9501_10_9_upl.5_800cells 

2.9 

PBX9501_10_9_upl.5_2000cells 

23.7 

8.2x 

PBX9501_10_34_upl.5_800cells 

2.6 

PBX9501  10  34  upl.5  2000cells 

17.0 

6.5x 

PBX9501  20  9  upl.5  800cells 

1.5 

PBX9501_20_9_upl.5_2000cells 

11.5 

7.7x 

The  average  percent  difference  for  all  the  detonation  properties,  except  tempera¬ 
ture,  is  1.5%.  The  average  data  output  for  2000  cell  simulations  is  7.5  times  larger 
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than  for  the  800  cell  simulations.  7.5  times  more  data  storage  is  not  worth  the  2% 
difference  in  the  detonation  properties  for  2000  computational  cells.  The  rest  of  the 
simulations  will  use  800  computational  cells  as  a  starting  grid  resolution  for  each 
metalized  case.  If  the  grid  resolution  does  not  provide  good  results  for  the  metalized 
case,  then  2000  computational  cells  will  be  utilized.  In  order  to  reduce  the  storage 
amount  the  data  saving  rate  can  be  changed.  The  data  saving  rate  is  defaulted  to 
save  a  fort  hie  every  10  time  steps.  For  2000  cells  the  rate  will  be  changed  to  save  a 
fort  hie  every  100  time  steps.  The  fewer  time  step  hies  saved  will  provide  less  data 
points  and  might  result  in  poor  results.  This  issue  will  be  considered  when  all  the 
metalized  explosive  simulation  data  at  800  computational  cells  is  studied  in  Chapter 
IV. 


3.6.2  PBX9501/AI  at  3  Different  Piston  Speeds 

Several  metalized  explosive  simulations  were  performed  at  three  diherent  piston 
speeds  to  test  the  limitations  of  the  MPEXS  code.  The  piston  speed  that  were 
considered  were  Up  =  0.5,  1.0,  and  1.5  km/s.  To  test  the  code  three  cases  were 
considered;  PBX9501_10_9,  PBX9501_10_13,  and  PBX9501_20_9.  This  covered 
simulations  at  the  lowest  concentration  of  aluminum,  which  is  10%  Al  mass  fraction, 
at  two  diherent  Al  particle  diameters.  It  also  covered  simulations  at  the  highest 
concentration  of  aluminum,  20%  Al  mass  fraction. 

The  piston  speed  of  Up  =  2.0  km/s  was  tested  for  the  three  cases.  The  simulation 
ran  with  no  issues  for  PBX9501_10_9  and  PBX9501_10_13  but  failed  immedi¬ 
ately  after  starting  for  PBX9501_20_9.  The  MPEXS  code  has  problems  performing 
calculations  with  high  aluminum  mass  fractions  and  high  piston  speeds.  A  reasons 
might  be  that  at  high  impact  speeds  the  explosive  mixture  might  have  compressibility 
issues.  This  high  metal  mass  fraction,  high  impact  speed  problem  has  been  reported 
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to  AFRL,  so  it  can  be  resolved  in  further  iterations  of  the  code. 

Since  PBX9501_20_9  failed  at  a  piston  speed  of  Up  =  2.0  km/s,  the  three  cases 
were  run  at  a  lower  piston  speed  of  Up  =  1.0  km/s.  The  MPEXS  code  ran  with  no 
issues  for  three  cases.  The  problem  with  this  lower  piston  speed  is  that  the  detonation 
wave  did  not  reach  steady  state  for  the  chosen  position  domain. 

Since  Up  =  2.0  km/s  was  too  fast  and  Up  =  1.0  km/s  was  too  slow,  a  piston  speed 
of  Up  =  1.5  km/s  was  tested  for  the  three  cases.  The  three  cases,  PBX9501_10_9, 
PBX9501_10_13,  and  PBX9501_20_9,  ran  with  no  issues  and  the  detonation  wave 
reached  steady  state  for  the  chosen  domain.  The  piston  speed  of  Up  =  1.5  km/s 
provided  good  results  and  will  be  utilized  for  all  the  simulations  what  will  be  discussed 
in  the  results  chapter. 

3.7  Design  of  Experiments 

Design  of  Experiments  (DOE)  is  the  process  of  determining  the  cause-and-effect 
relationships  between  input  parameters  and  how  these  relationships  affect  the  output. 
The  purpose  for  such  analysis  is  to  optimize  the  output  by  looking  at  the  inputs.  A 
stochastic  experiment  has  inputs  that  are  uncontrollable  and  are  usually  not  solved 
analytically.  DOE  methods  are  useful  for  stochastic  experiments  to  determine  the 
effect  that  these  uncontrollable  inputs  have  on  the  output.  In  deterministic  experi¬ 
ments  the  inputs  are  controllable  and  the  relationship  between  input  and  output  is 
conclusively  determined.  The  input  parameters  for  this  research  are  particle  size  and 
mass  fraction  and  both  are  controllable  and  are  independent  of  each  other.  The  six 
different  aluminum  particle  sizes  are  discrete  because  these  sizes  are  determined  by 
the  manufacturer.  The  aluminum  mass  fraction  is  determined  by  the  user.  These 
input  parameters  are  controllable  and  independent  of  each  other,  therefore  this  re¬ 
search  is  considered  deterministic.  DOE  methods  were  considered  for  the  multiphase 
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code  simulation  runs  but  there  was  no  benefit  in  using  a  DOE  method. 

3.8  Summary 

Section  3.2  explained  the  different  modeling  methods  that  lead  to  multiphase 
computational  modeling.  After  the  modeling  method  for  multiphase  explosives  was 
described,  a  brief  summary  was  given  of  the  main  hies  utilized  in  the  MPEXS  code  in 
Section  3.3.  The  MPEXS  code  will  be  utilized  to  perform  all  the  simulations  of  metal- 
ized  PBX9501  explosive.  Section  3.4  describes  the  different  simulations  runs  that  will 
be  performed  at  different  Al  mass  fractions  and  Al  particle  diameters.  Section  3.5  and 
Section  3.6  explained  the  different  preliminary  studies  that  were  performed  to  hnd 
the  limitations  of  the  code  and  selected  the  initial  parameters  for  all  the  simulation 
cases. 
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IV.  Data  Discussion 


4.1  Introduction 

This  chapter  will  present  the  results  from  the  three  main  studies  that  were  per¬ 
formed.  The  purpose  of  these  studies  was  to  evaluate  the  stability  and  limitations 
of  the  MPEXS  code.  The  first  study  is  a  short  study  with  neat  PBX9501  to  com¬ 
pare  run-to-detonation  properties  against  a  single  phase  explosive  code  developed  by 
AFRL.  This  study  will  verify  that  the  MPEXS  code  outputs  reasonable  data.  The 
second  study  is  a  convergence  test  of  detonation  properties  by  varying  the  number 
of  computational  cells  for  metalized  PBX9501.  The  goal  of  this  study  was  to  test 
the  stability  of  the  code  by  looking  for  a  convergence  in  the  data  as  the  number  of 
computational  cells  increased,  decreasing  the  grid  size.  The  third  study  is  to  study 
the  effects  of  aluminum  on  PBX9501  by  performing  several  simulations  at  a  range  of 
aluminum  mass  fractions  and  aluminum  particle  diameters.  The  results  of  the  three 
studies  are  presented  in  the  three  sections  below. 

4.2  Run-to-Detonation  Properties  for  Neat  PBX9501 

The  MPEXS  code  was  developed  to  run  simulations  for  multiphase  explosives  that 
contain  at  least  two  solid  components  like  an  explosive  and  a  metal.  For  this  research, 
the  chosen  explosive  is  PBX9501  and  the  chosen  non-reactive  metal  is  aluminum. 
The  MPEXS  code  also  has  the  capability  to  run  single  phase  explosive  simulations. 
Before  multiphase  simulations  are  considered,  the  MPEXS  code  must  hrst  be  tested 
with  single  phase  explosive  simulations. 

The  MPEXS  code  was  utilized  to  perform  a  short  Run2Det  study  of  neat  PBX9501 
at  several  piston  speeds  Up  =  0.5,  Up  =  1.0,  and  Up  =  1.5  km/s.  The  descrip¬ 
tion  “Run-to-detonation”  is  abbreviated  to  “Run2Det”  in  this  report.  The  MPEXS 
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Run2Det  results  will  be  compared  to  the  results  of  the  Single-Phase  Explosive  Simu¬ 
lation  (SPEXS)  code.  The  SPEXS  code  was  developed  and  utilized  by  AFRL’s  High 
Explosives  Research  and  Development  (HERD)  group  to  determine  Run2Det  prop¬ 
erties  for  single  phase  explosives.  The  goal  of  this  study  is  to  compare  the  MPEXS 
results  to  SPEXS  results.  The  equations  of  state  utilized  for  the  MPEXS  code  are 
different  from  the  SPEXS  code,  therefore  the  Run2Det  properties  should  be  different. 
Although  the  results  will  differ,  the  MPEXS  data  should  display  a  similar  trend  as 
the  SPEXS  data.  A  qualitative  comparison  will  verify  the  MPEXS  code  is  working 
properly. 

Two  important  parameters  in  explosive  studies  are  run-to-detonation  time  and 
run-to-detonation  distance.  The  Run2Det  time  and  distance  is  the  time  and  loca¬ 
tion  the  explosive  detonation  transitions  to  steady  state.  Detonation  properties  in 
the  transition  stage  provide  information  with  little  use  because  these  properties  are 
constantly  changing  and  unstable.  The  Run2Det  time  versus  distance  diagrams  pro¬ 
vide  an  illustration  of  the  transition  from  the  initial  shock  wave  to  the  detonation 
shock  wave.  These  time  versus  distance  diagrams,  often  call  “pop-plots” ,  provide  the 
Run2Det  time  and  distance  as  a  function  of  shock  inpnt  pressnre.  The  threshold  input 
pressure  is  set  by  the  user  in  the  postprocess. f90  file  described  in  Section  3.3. 1.4. 
Pop-plots  are  utilized  to  determine  the  transition  time  and  location  to  a  steady  state 
detonation. 

There  are  two  slopes  connected  by  a  smooth  transition  in  the  pop-plot.  The  slope 
closest  to  the  zero  position  represents  the  initial  shock  wave.  During  the  initial  shock 
wave  the  detonation  is  still  in  development.  The  second  slope  represents  the  steady 
state  detonation  shock  wave.  The  transition  from  one  slope  to  the  other  is  very 
gradual  and  difhcnlt  to  detect.  The  goal  is  to  determine  the  time  and  position  when 
the  detonation  reaches  steady  state.  The  time  and  position  can  be  approximated  by 
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finding  a  point  where  the  second  slope  begins,  showing  a  constant  rate  of  change. 
Figure  21  to  Figure  23  show  the  pop-plots  for  the  three  piston  speeds  tested  in  this 
study.  The  pop-plots  have  a  pressure  gradient  bar  on  the  right  side  of  the  plot.  The 
pressure  is  the  current  pressure  at  the  time  and  position  within  the  chart.  If  the  color 
gradient  that  represents  the  current  pressure  is  followed  along  the  slope  of  the  hrst 
slope,  it  is  seen  that  pressure  increases  as  the  initial  shock  wave  develops.  Along  the 
transition  to  the  steady  state  slope  line  the  pressure  increases  rapidly,  then  at  some 
point  the  pressure  stabilizes  to  a  constant  pressure.  The  point  where  the  pressure 
values  begin  to  stabilize  is  the  point  that  is  chosen  to  be  the  Run2Det  time  and 
position.  The  point  at  the  beginning  of  the  second  slope  can  be  seen  on  Figure  21  to 
Figure  23.  The  chosen  pressure  for  this  neat  PBX9501  study  is  3.0el0  Pa  (30  GPa). 

Three  pop-plots  were  generated  for  the  three  simulation  runs  at  piston  speeds  of 
Up  =  0.5,  Up  =  1.0,  and  Up  =  1.5  km/s.  Figure  21  is  the  pop-plot  for  Neat  PBX9501 
at  Up  =  0.5  km/s  and  shows  a  Run2Det  time  of  2.1705  ps  and  distance  of  5.001 
millimeters.  Figure  22  is  the  pop-plot  for  Neat  PBX9501  at  Up  =  1.0  km/s  and  shows 
a  Run2Det  time  of  0.41043  ps  and  distance  of  1.4125  millimeters.  Figure  23  is  the 
pop-plot  for  Neat  PBX9501  at  Up  =  1.5  km/s  and  shows  a  Run2Det  time  of  0.15460 
fis  and  distance  of  0.60063  millimeters. 
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Figure  21.  Neat  PBX9501  Pop-Plot  {up  —  0.5  L  =  18  Nx  =  3000) 
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Figure  22.  Neat  PBX9501  Pop-Plot  {up  —  1.0  L  =  4  Nx  =  800) 
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Figure  23.  Neat  PBX9501  Pop-Plot  {up  —  1.5  L  =  1  Nx  =  800) 

Table  8  is  a  summary  of  all  the  SPEXS  and  MPEXS  simulation  results  for  the 
three  piston  speeds.  As  seen  on  Table  8,  Run^^l  failed  for  the  simulations  run  with 
the  MPEXS  code.  The  reason  for  the  code  failure  was  that  the  grid  size  was  too 
large.  The  grid  size  is  calculated  by  taking  the  total  length  divided  by  the  number 
of  computational  cells  utilized.  The  length  for  Run^l  was  18  mm  and  the  number 
of  cells  utilized  was  800  cells  resulting  in  a  grid  size  of  0.022  mm/cell.  The  grid  size 
was  increased  by  using  3000  cells  over  the  same  length  for  Run7^2,  resulting  in  a  grid 
size  of  0.006  mm/ cell.  The  results  for  both  codes  are  as  follows. 
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Table  8.  Neat  PBX9501  Run-to-Detonation  Properties 


Grid 

SPEXS 

MPEXS 

SPEXS 

MPEXS 

Run  # 

Up  (km/s) 

Size 

Run2Det 

Run2Det 

Run2Det 

Run2Det 

{mm/cell) 

X  (mm) 

X  (mm) 

Time  (fis) 

Time  (fis) 

1 

0.5 

0.022 

15.00 

Failed 

4.50 

Failed 

2 

0.5 

0.006 

15.00 

5.00 

4.50 

2.17 

3 

1.0 

0.005 

2.00 

1.41 

0.45 

0.41 

4 

1.5 

0.001 

0.30 

0.60 

0.10 

0.15 

Figure  24  is  a  bar  chart  to  compare  the  Run2Det  distance  for  both  the  SPEXS 
and  MPEXS  codes.  The  Run2Det  distance  for  each  piston  speed  are  shown  in  blue 
for  SPEXS  and  red  for  MPEXS.  As  the  piston  speed  increases  the  Run2Det  distance 
decreases.  The  average  rate  of  change  between  Up  =  0.5  and  Up  =  1.0  is  26  mm/ {km/ s) 
for  SPEXS  and  7.18  mm/{km/s)  for  MPEXS.  The  average  rate  of  change  between 
Up  =  1.0  and  Up  =  1.5  is  3.40  mm/{km/s)  for  SPEXS  and  1.62  mm/{km/s)  for 
MPEXS.  The  reason  for  the  difference  in  results  between  both  codes  is  due  to  the  use 
of  different  equations  of  state  for  the  reaction.  The  average  rate  of  change  between 
each  set  of  points  is  larger  for  the  SPEXS  code  bnt  the  takeaway  is  that  the  MPEXS 
code  follows  the  same  trend. 
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Figure  24.  NeatPBX9501  R2D  Distance  Chart 


Figure  25  shows  the  Run2Det  time  at  different  piston  speeds  for  both  the  SPEXS 
and  MPEXS  codes.  As  the  piston  speed  increases  the  Run2Det  time  decreases.  The 
average  rate  of  change  between  Up  =  0.5  and  Up  =  1.0  is  8.10  ns/ikm/s)  for  SPEXS 
and  3.52  fis/{km/s)  for  MPEXS.  The  average  rate  of  change  between  Up  =  1.0  and 
Up  =  1.5  is  0.70  jjis/{km/s)  for  SPEXS  and  0.52  ns/{km/s)  for  MPEXS.  The  average 
rate  of  change  between  each  set  of  points  is  larger  for  the  SPEXS  code.  Again  the 
difference  in  the  results  between  both  codes  is  that  different  equations  of  state  were 
utilized. 
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Figure  25.  NeatPBX9501  R2D  Time  Chart 


4.3  PBX9501/AL  Convergence  Test 

Another  study  that  was  performed  was  a  convergence  test  of  detonation  properties 
for  different  computational  cell  sizes.  The  computational  cell  sizes  play  an  important 
role  in  the  grid  resolution  of  the  simulation.  The  number  of  computational  cells, 
Nx,  per  unit  length  is  the  grid  size.  The  MPEXS  code  was  designed  to  perform 
simulations  at  the  macroscopic  scale.  If  the  grid  size  is  large  the  code  will  stop 
because  the  equations  can  not  handle  macro  scale  calculations.  If  the  grid  size  is 
very  small  then  the  code  will  run  into  problems  were  the  grid  size  is  smaller  than 
the  particle  diameters.  The  convergence  test  was  performed  with  PBX9501_10_9 
at  a  piston  speed  of  1.0  km/s.  The  computational  cell  sizes  were  varied  from  200 
to  3000  cells.  There  are  three  detonation  properties  that  were  considered  for  this 
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convergence  test,  Run2Det  time,  Run2Det  distance,  and  detonation  velocity.  The 
Run2Det  distance  and  time  are  approximated  by  analyzing  the  pop-plot  chart  and 
hnding  the  point  were  the  detonation  transitions  to  steady-state  as  described  earlier. 
The  point  was  chosen  at  the  beginning  of  the  steady-state  slope  at  a  pressure  value 
of  35  GPa  for  all  metalized  explosive  cases. 

Table  9  summarizes  the  results  for  the  convergence  test  for  the  computational  cell 
range.  Run  1  and  2  failed  due  to  the  grid  size  being  too  large.  Based  on  this  study 
and  the  neat  explosive  study  from  Section  4.2  it  was  concluded  that  for  any  grid 
size  larger  than  9.9e-3  mm/cell  will  make  the  code  fail.  The  MPEXS  code  ran  for 
simulations  with  800  or  more  computational  cells.  The  Run2Det  distance  and  time, 
and  detonation  velocity  results  are  listed  below,  in  Table  9. 

The  hnal  stop  time  for  all  the  runs  in  this  study  was  set  at  thnal  =  2.0  /is  in  the 
maininput.txt  hie.  It  was  observed  that  all  the  simulations  did  not  run  to  the  hnal 
run  time.  The  runs  stopped  around  1.2  /xs  and  runs  10  and  11  stopped  around  1.4 
fis.  An  explanation  for  the  early  stop  time  is  that  the  errors  occurring  in  the  code 
are  due  to  larger  numerical  error  for  lower  resolutions  simulations.  As  the  resolution 
increases,  the  simulations  run  longer  due  to  a  decrease  in  numerical  error.  The  early 
stop  time  issue  will  be  reported  to  the  HERD  to  be  corrected  in  future  code  versions. 
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Table  9.  PBX9501  10  9  Convergence  Test 


Run  7^ 

Nx 

Grid  Size 

{mm/ cell) 

Run2Det 

X  (mm) 

Run2Det 

Time  (/is) 

Det. 

Vel.  (km/s) 

1 

200 

0.0300 

Failed 

Failed 

Failed 

2 

400 

0.0150 

Failed 

Failed 

Failed 

3 

800 

0.0075 

3.761 

1.193 

6.142 

4 

1000 

0.0060 

3.699 

1.182 

6.576 

5 

1300 

0.0046 

3.722 

1.182 

7.013 

6 

1600 

0.0037 

3.733 

1.183 

7.251 

7 

1800 

0.0033 

3.738 

1.184 

7.387 

8 

2000 

0.0030 

3.743 

1.184 

7.667 

9 

2500 

0.0024 

3.757 

1.186 

7.957 

10 

3000 

0.0020 

3.739 

1.183 

7.940 

The  results  for  the  Run2Det  distance  are  plotted  on  Figure  26.  For  smaller  com¬ 
putational  cells  the  Run2Det  distance  varies  but  as  the  cell  number  increases  the 
Run2Det  distance  seems  to  stabilize  at  a  distance  around  3.74  mm. 
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Figure  26.  Convergence  Test  for  PBX9501_10_9  Run2Det  Distance 


Figure  27  is  a  plot  of  Run2Det  time  in  microseconds  versus  computational  cells. 
The  Run2Det  time  has  the  same  relationship  to  Nx  as  Run2Det  distance.  At  smaller 
cell  amounts  the  results  are  unstable  but  as  the  number  of  cells  increases  the  Run2Det 
time  stabilizes  around  1.183  ps. 
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Figure  27.  Convergence  Test  for  PBX9501  10  9  Run2Det  Time 
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Figure  28  plots  the  detonation  velocity  in  kilometers  per  second  versus  compu¬ 
tational  cells.  As  the  resolution  increases  the  detonation  velocity  results  stabilize 
around  7.94  km/ s.  The  results  for  Run2Det  distance,  Run2Det  time,  and  detonation 
velocity  behave  as  expected  where  the  results  should  converge  to  a  steady  value  as 
the  resolution  is  increased. 


Detonation  Velocity 
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Figure  28.  Convergence  Test  for  PBX9501_10_9  Detonation  Velocity 


Figure  29  through  Figure  31  are  the  pop-plots  for  three  different  computational 
cells,  800,  1600,  and  3000  cells,  for  PBX9501_10_9.  The  domain  length  for  each 
simulation  is  from  the  left  boundary  location  (0  mm)  to  the  right  boundary  location 
(6  mm).  All  the  pop-plots  are  very  similar  since  Run2Det  distance  has  a  standard 
deviation  of  0.019  mm,  a  standard  deviation  of  0.003  fis  for  Run2Det  time,  and  a 
standard  deviation  of  0.644  km/s  for  detonation  velocity.  Figure  31  looks  different 
than  the  other  simulations  due  to  the  longer  running  time  when  3000  computational 
cells  were  utilized,  but  the  Run2Det  point  is  almost  the  same  for  the  three  plots. 
Although  the  simulations  stopped  earlier  than  the  hnal  simulation  time,  the  data 
collected  was  enough  to  determine  steady  state  detonation  properties  like  Run2Det 
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distance,  Run2Det  time,  and  detonation  velocities. 
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Figure  29.  PBX9501_10_9  Pop-Plot  (Nx=800  cells) 
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Figure  30.  PBX9501  10  9  Pop-Plot  (Nx=1600  cells) 
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Figure  31.  PBX9501  10  9  Pop-Plot  (Nx=3000  cells) 

The  convergence  test  was  performed  to  determine  the  stability  of  grid  convergence 
and  the  limits  of  the  MPEXS  code.  The  grid  convergence  of  the  code  was  proved  by 
the  converging  results  as  the  number  of  computational  cells  increased.  The  number 
of  computational  cells  is  related  to  the  resolution  or  the  grid  size  of  the  simulation. 
The  hner  the  resolution  the  more  precise  the  results  and  should  converge  as  seen  in 
this  study.  The  limits  of  the  MPEXS  were  also  determined  during  this  study.  The 
MPEXS  code  did  not  run  with  computational  sizes  such  as  Nx  =  200  and  Nx  =  400. 
For  Nx  =  800  and  Nx  =  1000  the  results  were  unpredictable. 

4.4  All  Cases  of  PBX9501/A1  to  Steady  State  Detonation 

The  third  study  was  performed  with  all  the  cases  of  PBX9501/A1,  which  covered  a 
wide  range  of  aluminum  particle  diameter  and  aluminum  mass  fraction.  The  six  alu¬ 
minum  diameters  were  tested  for  each  of  the  three  aluminum  mass  fractions,  resulting 
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in  18  simulation  runs.  For  example,  the  aluminum  mass  fraction  of  10  percent  was 
tested  with  9,  13,  19,  23,  29,  and  34  micron  diameter.  The  same  six  diameters  were 
utilized  for  15  and  20  percent  aluminum  mass  fraction.  The  18  runs  were  performed 
for  two  computational  cell,  800  and  2000,  for  a  total  of  36  simulation  runs.  These 
two  computational  cell  sizes  of  Nx  =  800  and  Nx  =  2000  were  chosen  to  compare 
the  data  for  both  a  low  and  high  grid  size.  A  piston  speed  of  Up  =  1.5  km/s  was  se¬ 
lected  for  all  the  PBX9501/A1  cases.  The  main  parameters  considered  were  Run2Det 
distance,  Run2Det  time,  detonation  velocity,  pressure,  temperature,  particle  velocity, 
and  volume  fraction. 

This  study  will  evaluate  the  stability  and  limits  of  the  MPEXS  code.  The  results 
will  be  analyzed  and  qualitative  conclusions  will  be  made  about  the  effects  of  Al  on 
PBX9501.  As  discussed  in  earlier  chapters,  comparing  simulation  data  to  empirical 
data  does  not  provide  any  benefit  to  this  research  due  to  the  complexity  of  plastically 
bonded  explosives  like  PBX9501. 

4.4.1  Run-to-Detonation  Properties  and  Detonation  Velocity 

The  18  runs  were  performed  at  800  computational  cells  at  Up  =  1.5  km/s.  The 
results  for  Run2Det  distance  and  time,  and  detonation  velocity  are  presented  on 
Table  10.  As  the  Al  mass  fraction  increases  from  10  to  15  to  20  percent,  the  Run2Det 
distance  decreases.  As  the  Al  diameter  increases  from  9  to  34  fim,  the  Run2Det 
distance  slightly  decreases.  The  effects  of  varying  the  mass  fraction  are  larger  than  the 
effects  of  varying  the  diameter  of  the  Al  for  Run2Det  distance.  As  the  mass  fraction 
is  increased  from  10  to  15  to  20  percent,  the  Run2Det  time  decreases  signihcantly. 
As  the  Al  diameter  is  increased  from  9  to  34  /rm  for  each  mass  fraction,  the  Run2Det 
time  varies  but  on  average  it  remains  constant.  The  results  for  detonation  velocity 
are  unpredictable  for  the  15  and  20  percent  cases. 
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The  same  18  runs  were  performed  again  for  2000  computational  cells  at  the  same 
piston  speed  Up  =  1.5  km/s.  The  results  for  Run2Det  distance  and  time,  and  deto¬ 
nation  velocity  are  presented  on  Table  11.  The  results  for  2000  computational  cells 
are  similar  from  the  800  computational  cells  cases  discussed  in  the  previous  para¬ 
graph  but  the  results  are  more  stable  for  detonation  velocity.  As  the  Al  mass  fraction 
increases  from  10  to  15  to  20  percent,  the  Run2Det  distance  decreases.  As  the  Al 
diameter  increases  from  9  to  34  /rm,  the  Run2Det  distance  slightly  decreases.  The 
effects  of  varying  the  mass  fraction  are  larger  than  the  effects  of  varying  the  diameter 
of  the  Al  for  Run2Det  distance.  As  the  mass  fraction  is  increased  from  10  to  15  to  20 
percent,  the  Run2Det  time  decreases  signihcantly.  As  the  Al  diameter  is  increased 
from  9  to  34  yum  for  each  mass  fraction,  the  Run2Det  time  varies  but  on  average  it 
remains  constant.  As  the  Al  mass  fraction  increases  from  10  to  15  to  20  percent,  the 
detonation  velocity  decreases.  Also,  as  the  Al  particle  diameter  increases  from  9  to 
34  yum  for  each  mass  fraction,  the  data  shows  the  detonation  velocity  decreases. 
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Table  10.  PBX9501/A1  Detonation  Properties  at  Up  —  1.5  kmjs  Nx  =  SOOcells 


Run  # 

Name 

Run2Det 

X  (mm) 

Run2Det 

Time  (/is) 

Det. 

Vel.  (km/s) 

1 

PBX9501_10_9 

1.016 

0.258 

8.024 

2 

PBX9501_10_13 

1.021 

0.259 

7.997 

3 

PBX9501_10_19 

1.026 

0.260 

7.978 

4 

PBX9501_10_23 

1.026 

0.261 

7.967 

5 

PBX9501_10_29 

1.021 

0.260 

7.954 

6 

PBX9501_10_34 

1.011 

0.258 

7.953 

7 

PBX9501_15_9 

0.939 

0.233 

7.749 

8 

PBX9501_15_13 

0.941 

0.234 

7.745 

9 

PBX9501_15_19 

0.939 

0.234 

7.621 

10 

PBX9501_15_23 

0.939 

0.235 

7.778 

11 

PBX9501_15_29 

0.929 

0.234 

7.807 

12 

PBX9501_15_34 

0.921 

0.233 

7.770 

13 

PBX9501_20_9 

0.866 

0.210 

7.378 

14 

PBX9501_20_13 

0.864 

0.211 

7.455 

15 

PBX9501_20_19 

0.869 

0.211 

7.490 

16 

PBX9501_20_23 

0.861 

0.212 

7.478 

17 

PBX9501_20_29 

0.851 

0.211 

7.460 

18 

PBX9501_20_34 

0.854 

0.211 

7.561 
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Table  11.  PBX9501/A1  Detonation  Properties  at  Up  —  1.5  kmjs  Nx  =  2000cells 


Run  # 

Name 

Run2Det 

X  (mm) 

Run2Det 

Time  (/is) 

Det. 

Vel.  (km/s) 

1 

PBX9501_10_9 

1.014 

0.255 

8.060 

2 

PBX9501_10_13 

1.020 

0.257 

8.026 

3 

PBX9501_10_19 

1.013 

0.257 

8.014 

4 

PBX9501_10_23 

1.012 

0.257 

8.005 

5 

PBX9501_10_29 

1.011 

0.257 

7.996 

6 

PBX9501_10_34 

1.017 

0.258 

7.997 

7 

PBX9501_15_9 

0.944 

0.232 

7.897 

8 

PBX9501_15_13 

0.943 

0.233 

7.866 

9 

PBX9501_15_19 

0.939 

0.233 

7.858 

10 

PBX9501_15_23 

0.938 

0.233 

7.852 

11 

PBX9501_15_29 

0.932 

0.233 

7.844 

12 

PBX9501_15_34 

0.926 

0.232 

7.848 

13 

PBX9501_20_9 

0.888 

0.214 

7.813 

14 

PBX9501_20_13 

0.887 

0.215 

7.815 

15 

PBX9501_20_19 

0.871 

0.212 

7.783 

16 

PBX9501_20_23 

0.873 

0.213 

7.765 

17 

PBX9501_20_29 

0.859 

0.211 

7.739 

18 

PBX9501_20_34 

0.850 

0.209 

7.744 
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Figure  32  and  Figure  33  shows  the  pop-plots  for  PBX9501  at  10%  aluminum  mass 
fraction  for  9  mm  and  34  mm  aluminum  diameters.  The  smallest  and  largest  of  Al 
particle  diameter  were  chosen  for  each  Al  mass  fraction.  The  pop-plots  for  the  Al 
particle  diameters  13,  19,  23,  29  are  very  similar  to  each  other  at  each  Al  mass  frac¬ 
tion.  For  brevity  of  the  report  all  the  pop-plots  were  not  included  but  the  results  are 
presented  on  Table  10  and  Table  11.  The  pop-plots  have  a  pressure  gradient  bar  on 
the  right  side  of  the  plot  whose  units  are  in  gigapascals.  The  Run2Det  properties  for 
PBX9501_10_9_upl.5_2000cells  are  Run2Det  time  =  0.255  ps  and  Run2Det  dis¬ 
tance  =  1.014  mm.  The  Run2Det  properties  for  PBX9501_10_34_upl.5_2000cells 
are  Run2Det  time  =  0.258  ps  and  Run2Det  distance  =  1.017  mm. 


0.5  1  1.5 

Position  (mm) 

Figure  32.  PBX9501_10_9  Pop-Plot  (Nx=2000  cells) 
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Figure  33.  PBX9501  10  34  Pop-Plot  (Nx=2000  cells) 

Figure  34  and  Figure  35  show  the  pop-plots  for  PBX9501  at  15%  aluminum  mass 
fraction  for  9  mm  and  34  mm  aluminum  diameters.  The  Run2Det  properties  for 
PBX9501_15_9_upl.5_2000cells  are  Run2Det  time  =  0.232  /rs  and  Run2Det  dis¬ 
tance  =  0.944  mm.  The  Run2Det  properties  for  PBX9501_15_34_upl.5_2000cells 
are  Run2Det  time  =  0.232  /is  and  Run2Det  distance  =  0.926  mm. 
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Figure  34.  PBX9501  15  9  Pop-Plot  (Nx=2000  cells) 
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Figure  35.  PBX9501_15_34  Pop-Plot  (Nx=2000  cells) 


Figure  36  and  Figure  37  show  the  pop-plots  for  PBX9501  at  20%  aluminum  mass 
fraction  for  9  fim  and  34  pm  aluminum  diameters.  Figure  36  seems  to  have  a  greater 
slope  than  the  rest  of  the  pop-plots  but  it  only  looks  different  due  to  the  shorter  run 
time  for  the  PBX9501_20_9  simulation  run.  The  results  for  Run2Det  distance  and 
time  for  both  cases  are  within  0.012  mm  and  0.001  ps  of  each  other.  The  Run2Det 
properties  for  PBX9501_20_9_upl.5_2000cells  are  Run2Det  time  =  0.214  ps  and 
Run2Det  distance  =  0.888  mm.  The  Run2Det  properties  for 

PBX9501_20_34_upl.5_2000cells  are  Run2Det  time  =  0.209  ps  and  Run2Det  dis¬ 
tance  =  0.850  mm. 


Position  (mm) 

Figure  36.  PBX9501  20  9  Pop-Plot  (Nx=2000  cells) 
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Figure  37.  PBX9501  20  34  Pop-Plot  (Nx=2000  cells) 


Figure  38  shows  the  Run2Det  distance  for  all  the  runs  at  both  computational  cell 
amounts,  800  and  2000  cells.  Runs  1-6  are  for  PBX9501  at  10%  Al  mass  fraction  for 
the  six  Al  particle  diameters.  The  red  square  data  points  are  for  Nx2000  and  the  blue 
diamond  data  points  are  for  Nx800.  The  other  runs  for  15  and  20  percent  are  also 
plotted  and  can  be  identihed  using  the  key  on  the  right  side  of  the  chart  in  Figure  38. 

The  Run2Det  distance  for  10%  PBX9501  at  Nx2000,  represented  by  the  red 
squares,  goes  up  and  down,  and  has  no  noticeable  trend.  The  Run2Det  distance 
for  10%  PBX9501  at  Nx800,  represented  by  the  blue  diamonds,  has  a  small  convex 
trend  where  the  distance  increases  and  then  decreases  as  the  Al  particle  diameter  in¬ 
creases.  The  Run2Det  distance  for  15%  PBX9501  at  Nx2000  and  Nx800,  represented 
by  the  purple  X  and  green  triangle  respectively,  both  have  a  slight  decrease  in  dis¬ 
tance  as  the  Al  particle  diameter  increases.  The  Run2Det  distance  for  20%  PBX9501 
at  Nx2000  and  Nx800,  represented  by  the  orange  circle  and  teal  star  respectively, 
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both  have  steady  values  then  decrease  slightly  and  then  steady  again.  In  general  a 
conclusion  can  be  made  that  as  the  aluminum  mass  fraction  increases  from  10  to  15 
to  20  percent,  the  Run2Det  distance  decreases.  As  the  aluminum  diameter  increases 
from  9  to  34  pm,  the  Run2Det  distance  slightly  decreases.  The  effects  of  varying  the 
mass  fraction  are  larger  than  the  effects  of  varying  the  diameter  of  the  aluminum. 


Run-to-Detonation  Distance 
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Figure  38.  PBX9501  All  Cases  Run2Det  Distance 


Figure  39  shows  the  Run2Det  time  for  all  the  runs  at  both  computational  cell 
amounts,  800  and  2000  cells.  The  Run2Det  time  for  10%  PBX9501  at  Nx2000, 
represented  by  the  red  squares,  increases,  decreases,  and  then  increases  again  not 
showing  any  trend  in  the  data.  The  Run2Det  time  for  10%  PBX9501  at  Nx800, 
represented  by  the  blue  diamonds,  has  a  small  convex  trend  as  the  aluminum  size 
increases.  The  Run2Det  time  for  15%  PBX9501  at  Nx2000  and  Nx800,  represented 
by  the  purple  X  and  green  triangle  respectively,  both  have  a  small  convex  trend  as 
the  Al  particle  diameter  increases.  The  results  for  both  case  are  almost  identical 
and  points  are  almost  on  top  of  each  other.  The  Run2Det  time  for  20%  PBX9501 


at  Nx2000,  represented  by  the  orange  circle,  has  a  slight  decrease  as  the  Al  particle 
diameter  increases.  The  Rnn2Det  time  for  20%  PBX9501  at  Nx800,  represented  by 
the  teal  star,  is  steady  as  the  Al  particle  diameter  increases.  The  effects  of  varying 
the  mass  fraction  are  larger  than  the  effects  of  varying  the  diameter  of  the  aluminum. 
When  the  mass  fraction  is  increased  from  10  to  15  to  20  percent,  the  Run2Det  time 
decreases  signihcantly.  When  the  aluminum  diameter  is  increased  from  9  to  34  pm 
for  each  mass  fraction,  the  Run2Det  time  varies  but  the  variation  is  small. 
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Figure  39.  PBX9501  All  Cases  Run2Det  Time 


Figure  40  plots  the  detonation  velocity  versus  the  run  numbers  at  the  three  differ¬ 
ent  aluminum  mass  fractions.  The  detonation  velocity  for  10%  PBX9501  at  Nx2000 
and  Nx800,  represented  by  the  red  squares  and  blue  diamonds  respectively,  both  have 
a  steady  decrease  as  the  Al  particle  diameter  increases.  The  detonation  velocity  for 
15%  PBX9501  at  Nx2000,  represented  by  the  purple  X,  has  a  steady  decrease  as  the 
Al  particle  diameter  increases.  The  detonation  velocity  for  15%  PBX9501  at  Nx800, 
represented  by  the  green  triangle,  has  a  very  unstable  pattern  with  large  deviations 
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in  the  detonation  velocity  as  the  Al  size  increases.  The  detonation  velocity  for  20% 
PBX9501  at  Nx2000,  represented  by  the  orange  circle,  has  a  steady  decrease  as  the 
Al  particle  diameter  increases.  The  detonation  velocity  for  20%  PBX9501  at  Nx800, 
represented  by  the  teal  star,  has  a  unstable  pattern  in  the  data.  For  this  case,  the 
detonation  velocity  increases  as  the  particle  diameter  increases,  which  is  opposite  of 
what  is  expected.  Considering  only  the  more  stable  results  using  2000  cells  a  few 
conclusions  can  be  made  for  detonation  velocity.  As  the  Al  mass  fraction  increases 
from  10  to  15  to  20  percent,  the  detonation  velocity  decreases.  Also,  as  the  Al  par¬ 
ticle  diameter  increases  from  9  to  34  /rm  for  each  mass  fraction,  the  data  shows  the 
detonation  velocity  decreases. 
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Figure  40.  PBX9501AllCases  Detonation  Velocity 


It  was  concluded  that  the  simulations  using  2000  computational  cells  provide  more 
consistent  results  for  both  varying  cases,  Al  mass  fraction  and  Al  particle  diameter. 
The  purpose  of  running  simulations  with  both  sizes  of  computational  cells,  800  and 
2000  cells,  is  study  the  effects  of  grid  size  on  metalized  explosive  simulations.  In 
Section  4.3,  it  was  found  that  at  lower  resolutions  the  detonation  property  results  are 


91 


unstable.  The  same  stability  issues  at  low  resolution  was  also  seen  for  the  metalized 
explosive  simulations  presented  in  this  section. 

4.4.2  Other  Detonation  Properties 

Detonation  properties  like  Run2Det  distance  and  time,  and  detonation  velocity  are 
provide  more  useful  information  about  explosion  but  other  properties  like  pressure, 
particle  velocity,  volume  fraction,  and  temperature  must  not  be  neglected.  Pressure 
of  the  shock  wave  provides  information  about  the  strength  of  the  shock  wave.  Parti¬ 
cle  velocity  provides  information  on  the  velocity  of  the  shock  wave.  Volume  fraction 
provides  information  about  the  chemical  reaction.  Temperature  provides  information 
about  the  chemical  reaction  but  as  described  earlier,  the  computational  model  does 
not  accurately  predict  the  actual  temperatures  near  the  shock  wave  due  to  the  insuf¬ 
ficient  resolution  near  the  reaction  zone,  therefore  those  plots  will  not  be  presented 
[5]. 

For  Figure  41  to  Figure  49,  pressure,  particle  velocity,  and  volume  fraction  will  be 
plotted  against  position.  These  hgures  are  a  snap-shoot  of  the  shock  wave  at  a  time 
after  the  wave  has  reached  steady-state.  To  compare  the  different  cases  to  each  other 
a  position  of  1.6  mm  was  chosen  for  all  the  plot.  Although  the  time  might  be  different 
for  each  case,  the  chosen  position  is  far  enough  that  steady-state  has  been  reached 
for  all  the  cases.  Once  steady-state  is  reached  all  the  properties  will  not  change  for 
the  range  being  considered  here. 

Figure  41  to  Figure  43  plot  pressure  versus  position  for  the  smallest  and  largest 
Al  particle  diameter,  9  and  34  pm,  at  each  Al  mass  fraction.  The  most  important 
pressure  for  these  plots  is  the  pressure  of  gas  reaction  product.  The  pressure  of  the 
metal  particles  is  not  important  because  it  is  not  driving  the  shock  wave.  The  pressure 
of  the  explosive  is  the  same  as  the  gas  pressure  until  the  explosive  volume  fraction 
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goes  to  zero.  The  data  for  the  pressure  of  the  explosive  is  irrelevant  after  the  volume 
fraction  goes  to  zero  because  all  the  explosive  material  has  been  consumed  by  the 
reaction. 

Figure  41a  shows  the  peak  gas  pressure  of  the  shock  wave  at  39.20  GPa  for 
PBX9501_10_9_upl.5.  Figure  41b  shows  the  peak  gas  pressure  of  the  shock  wave 
at  43.64  GPa  for  PBX9501_10_34_upl.5.  It  is  observed  that  at  10%  Al  mass 
fraction  the  shock  wave  pressure  increases  as  the  Al  particle  diameter  increases. 


Pressure  vs  Position 


Pressure  vs  Position 


(a)  PBX9501_10_9_upl.5_2000cells 


(b)  PBX9501_10_34_upl.5_2000cells 


Figure  41.  PBX9501  10%  Al  Mass  Fraction  Pressure  Plot 


Figure  42a  and  Figure  42b  plot  pressure  versus  position  for  PBX9501_15_9_upl.5 
and  PBX9501_15_34_upl.5.  Figure  42a  shows  the  peak  gas  pressure  of  the  shock 
wave  at  36.51  GPa  for  PBX9501_15_9_upl.5.  Figure  42b  shows  the  peak  gas  pres¬ 
sure  of  the  shock  wave  at  42.31  GPa  for  PBX9501_15_34_upl.5.  It  is  observed 
that  at  15%  Al  mass  fraction  the  shock  wave  pressure  increases  as  the  Al  particle 
diameter  increases. 
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Pressure  vs  Position 


Pressure  vs  Position 


(a)  PBX9501_15_9_upl.5_2000cells  (b)  PBX9501_15_34_upl.5_2000cells 

Figure  42.  PBX9501  15%  Al  Mass  Fraction  Pressure  Plot 

Figure  43a  and  Figure  43b  plot  pressure  versus  position  for  PBX9501_20_9_upl.5 
and  PBX9501_20_34_upl.5.  Figure  43a  shows  the  location  of  shock  wave  to  be  at 
around  1.15  mm.  The  reason  for  the  shock  wave  stopping  at  this  location  is  due  to 
the  simulation  quiting  at  t  =  0.257  ns.  Even  though  the  simulation  stopped  early, 
the  detonation  reached  steady-state  prior  to  the  quiting  time  and  the  results  are  still 
valid.  Figure  43a  shows  the  peak  gas  pressure  of  the  shock  wave  at  31.23  GPa  for 
PBX9501_20_9_upl.5.  Figure  43b  shows  the  peak  gas  pressure  of  the  shock  wave  at 
37.61  GPa  for  PBX9501_20_34_upl.5.  It  is  observed  that  at  20%  Al  mass  fraction 
the  shock  wave  pressure  increases  as  the  Al  particle  diameter  increases. 
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Pressure  vs  Position 


Pressure  vs  Position 


(a)  PBX9501_20_9_upl.5_2000cells  (b)  PBX9501_20_34_upl.5_2000cells 

Figure  43.  PBX9501  20%  Al  Mass  Fraction  Pressure  Plot 

Figure  44  to  Figure  46  plot  particle  velocity  versus  position  for  the  smallest  and 
largest  Al  particle  diameter,  9  and  34  /rm,  at  each  Al  mass  fraction.  The  most 
important  particle  velocity  for  these  plots  is  the  particle  velocity  of  gas  reaction 
product.  The  data  for  the  particle  velocity  of  the  explosive  is  irrelevant  after  the 
volume  fraction  goes  to  zero  because  all  the  explosive  material  has  been  consumed 
by  the  reaction.  Figure  44a  shows  the  peak  gas  particle  velocity  of  the  shock  wave 
at  3.609  km/s  for  PBX9501_10_9_upl.5.  Figure  44b  shows  the  peak  gas  particle 
velocity  of  the  shock  wave  at  3.727  km/s  for  PBX9501_10_34_upl.5.  It  is  observed 
that  at  10%  Al  mass  fraction  the  shock  wave  particle  velocity  increases  slightly  as 
the  Al  particle  diameter  increases. 
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Particle  Velocity  vs  Position  Particle  Velocity  vs  Position 


(a)  PBX9501_10_9_upl.5_2000cells  (b)  PBX9501_10_34_upl.5_2000cells 

Figure  44.  PBX9501  10%  Al  Mass  Fraction  Velocity  Plot 

Figure  45a  and  Figure  45b  plot  particle  velocity  versus  position  for 
PBX9501_15_9_upl.5  and  PBX9501_15_34_upl.5.  Figure  45a  shows  the  peak 
gas  particle  velocity  of  the  shock  wave  at  3.580  km/s  for  PBX9501_15_9_upl.5. 
Figure  45b  shows  the  peak  gas  particle  velocity  of  the  shock  wave  at  3.774  km/s  for 
PBX9501_15_34_upl.5.  It  is  observed  that  at  15%  Al  mass  fraction  the  shock  wave 
particle  velocity  increases  slightly  as  the  Al  particle  diameter  increases. 
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Particle  Velocity  vs  Position  Particle  Velocity  vs  Position 


(a)  PBX9501_15_9_upl.5_2000cells  (b)  PBX9501_15_34_upl.5_2000cells 

Figure  45.  PBX9501  15%  Al  Mass  Fraction  Velocity  Plot 

Figure  46a  and  Figure  46b  plot  particle  velocity  versus  position  for 
PBX9501_20_9_upl.5  and  PBX9501_20_34_upl.5.  Figure  46a  shows  the  peak 
gas  particle  velocity  of  the  shock  wave  at  3.294  km/s  for  PBX9501_20_9_upl.5. 
Figure  46b  shows  the  peak  gas  particle  velocity  of  the  shock  wave  at  3.591  km/s  for 
PBX9501_20_34_upl.5.  It  is  observed  that  at  20%  Al  mass  fraction  the  shock  wave 
particle  velocity  increases  slightly  as  the  Al  particle  diameter  increases. 
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Particle  Velocity  vs  Position  Particle  Velocity  vs  Position 


(a)  PBX9501_20_9_upl.5_2000cells  (b)  PBX9501_20_34_upl.5_2000cells 

Figure  46.  PBX9501  20%  Al  Mass  Fraction  Velocity  Plot 


Figure  47  to  Figure  49  plot  volume  fraction  versus  position  for  the  smallest  and 
largest  Al  particle  diameter,  9  and  34  /rm,  at  each  Al  mass  fraction.  All  the  volume 
fraction  plots  presented  in  the  hgures  below  look  the  same.  The  reason  all  the  plots 
look  the  same  is  because  all  the  simulations  should  be  at  steady-state  at  this  position. 
The  volume  fraction  for  the  explosive  material,  represented  by  the  blue  line,  should 
go  to  zero  shortly  after  the  shock  wave  has  passed.  The  volume  fraction  for  the  gas, 
represented  by  the  red  line,  should  go  to  total  volume  fraction  set  by  the  user  after  the 
shock  wave  has  passed.  The  volume  fraction  for  the  metal,  represented  by  the  green 
line,  should  remain  constant  because  the  metal  in  non  reactive  in  the  simulation.  The 
explosive  volume  fraction  should  exponentially  decay  to  zero  very  sharply  and  the  gas 
volume  fraction  should  exponentially  increase  very  sharply  to  the  total  solids  volume 
fraction.  All  these  conditions  are  met  by  all  the  simulation  at  the  plotted  position 
and  time,  verifying  that  all  simulations  are  at  steady-state. 
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Volume  Fraction  Volume  Fraction 


1 


Volume  Fraction  vs  Position 


Volume  Fraction  vs  Position 


(a)  PBX9501_10_9_upl.5_2000cells  (b)  PBX9501_10_34_upl.5_2000cells 


Figure  47.  PBX9501  10%  Al  Mass  Fraction  Volume  Fraction  Plot 


Volume  Fraction  vs  Position  Volume  Fraction  vs  Position 


(a)  PBX9501_15_9_upl.5_2000cells 


(b)  PBX9501_15_34_upl.5_2000cells 


Figure  48.  PBX9501  15%  Al  Mass  Fraction  Volume  Fraction  Plot 
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Volume  Fraction  vs  Position  Volume  Fraction  vs  Position 


(a)  PBX9501_20_9_upl.5_2000cells 


(b)  PBX9501_20_34_upl.5_2000cells 


Figure  49.  PBX9501  20%  Al  Mass  Fraction  Volume  Fraction  Plot 


Figure  50a  shows  the  peak  gas  density  of  the  shock  wave  at  2319  kg/rn^  for 
PBX9501_10_9_upl.5.  Figure  50b  shows  the  peak  gas  density  of  the  shock  wave 
at  2419  kg/m^  for  PBX9501_10_34_upl.5.  It  is  observed  that  at  10%  Al  mass 
fraction  the  shock  wave  density  increases  as  the  Al  particle  diameter  increases. 


Density  vs  Position  Density  vs  Position 


(a)  PBX9501_10_9_upl.5_2000cells  (b)  PBX9501_10_34_upl.5_2000cells 

Figure  50.  PBX9501  10%  Al  Mass  Fraction  Density  Plot 


Figure  51a  and  Figure  51b  plot  density  versus  position  for  PBX9501_15_9_upl.5 
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and  PBX9501_15_34_upl.5.  Figure  51a  shows  the  peak  gas  density  of  the  shock 
wave  at  2243  kg/rrr’  for  PBX9501_15_9_upl.5.  Figure  51b  shows  the  peak  gas 
density  of  the  shock  wave  at  2361  kg/m?  for  PBX9501_15_34_upl.5.  It  is  observed 
that  at  15%  Al  mass  fraction  the  shock  wave  density  increases  as  the  Al  particle 
diameter  increases. 


Density  vs  Position 


Density  vs  Position 


(a)  PBX9501_15_9_upl.5_2000cells 


(b)  PBX9501_15_34_upl.5_2000cells 


Figure  51.  PBX9501  15%  Al  Mass  Fraction  Density  Plot 


Figure  52a  and  Figure  52b  plot  density  versus  position  for  PBX9501_20_9_upl.5 
and  PBX9501_20_34_upl.5.  Figure  52a  shows  the  peak  gas  density  of  the  shock 
wave  at  2168  kg/m^  for  PBX9501_20_9_upl.5.  Figure  52b  shows  the  peak  gas 
density  of  the  shock  wave  at  2270  kg/m^  for  PBX9501_20_34_upl.5.  It  is  observed 
that  at  20%  Al  mass  fraction  the  shock  wave  density  increases  as  the  Al  particle 
diameter  increases. 
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Density  vs  Position  Density  vs  Position 


(a)  PBX9501_20_9_upl.5_2000cells 


(b)  PBX9501_20_34_upl.5_2000cells 


Figure  52.  PBX9501  20%  Al  Mass  Fraction  Density  Plot 


Table  12  provides  the  data  collected  for  the  shock  wave  gas  pressure,  gas  products 
particle  velocity,  and  gas  density.  Run  ^  1  through  6  are  the  simulations  for  10%  Al 
mass  fraction  from  9  /tm  to  34  /rm  Al  particle  diameter.  Run  ^  7  through  12  are  the 
simulations  for  15%  Al  mass  fraction  from  9  pm  to  34  pm  Al  particle  diameter.  Run 
13  through  18  are  the  simulations  for  20%  Al  mass  fraction  from  9  pm  to  34  pm 
Al  particle  diameter.  The  gas  pressure  is  in  gigapascals,  the  gas  particle  velocity  is 
in  kilometers  per  second,  and  the  density  is  in  kilograms  per  meter  cubed. 
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Table  12.  PBX9501/A1  Gas  Properties  at  Up  —  1.5  km/s  Nx  =  2000cells 


Run 

# 

Name 

Pressure 

(GPa) 

Part.  Vel. 

{km/s) 

Density 

{kg/m?) 

1 

PBX9501_10_9 

39.200 

3.609 

2319 

2 

PBX9501_10_13 

40.640 

3.779 

2336 

3 

PBX9501_10_19 

41.300 

3.718 

2363 

4 

PBX9501_10_23 

41.960 

3.719 

2381 

5 

PBX9501_10_29 

42.480 

3.798 

2403 

6 

PBX9501_10_34 

43.640 

3.727 

2419 

7 

PBX9501_15_9 

36.510 

3.580 

2243 

8 

PBX9501_15_13 

37.700 

3.601 

2259 

9 

PBX9501_15_19 

38.940 

3.831 

2297 

10 

PBX9501_15_23 

40.530 

3.737 

2319 

11 

PBX9501_15_29 

41.330 

3.700 

2346 

12 

PBX9501_15_34 

42.310 

3.774 

2361 

13 

PBX9501_20_9 

31.230 

3.294 

2168 

14 

PBX9501_20_13 

33.600 

3.411 

2167 

15 

PBX9501_20_19 

34.680 

3.573 

2197 

16 

PBX9501_20_23 

35.850 

3.464 

2218 

17 

PBX9501_20_29 

36.600 

3.652 

2245 

18 

PBX9501_20_34 

37.610 

3.591 

2270 

The  results  from  the  table  above  are  presented  in  a  chart  to  better  observe  the 
trends  in  the  data.  Figure  53  plots  the  pressure  versus  Run  for  the  three  aluminum 
mass  fractions.  The  first  six  data  point  are  for  the  six  Al  particle  diameters  at  10%  Al 
mass  fraction  and  are  represented  by  the  blue  diamonds.  The  next  six  data  point  are 
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for  the  six  Al  particle  diameters  at  15%  Al  mass  fraction  and  are  represented  by  the 
red  squares.  The  last  six  data  point  are  for  the  six  Al  particle  diameters  at  20%  Al 
mass  fraction  and  are  represented  by  the  green  triangles.  It  can  be  seen  from  Figure 
53  the  gas  pressure  steadily  increases  as  the  Al  particle  diameter  increases  at  10%  Al 
mass  fraction.  For  15%  Al  mass  fraction  the  hrst  pressure  is  less  than  than  all  the 
10%  Al  mass  fraction  cases.  The  same  trend  is  observed  for  15%  Al  mass  fraction, 
where  the  gas  pressure  steadily  increases  as  the  Al  particle  diameter  increases.  Again 
for  20%  Al  mass  fraction  the  hrst  pressure  is  lower  than  all  the  10%  and  15%  Al  mass 
fraction  cases.  The  same  trend  is  observed  for  20%  Al  mass  fraction  cases  where  the 
gas  pressure  steadily  increases  as  the  Al  particle  diameter  increases. 


PBX9501/A1  Pressure 
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Figure  53.  PBX9501  All  Cases  Gas  Pressure 


Figure  54  plots  the  gas  particle  velocity  versus  Run  ^  for  the  three  aluminum  mass 
fractions.  The  same  color  and  symbols  are  utilized  for  the  different  simulation  runs 
as  in  Figure  53.  All  the  18  simulations  have  the  same  gas  particle  velocity  pattern. 
As  the  Al  particle  diameter  increases  in  each  mass  fraction  case,  the  particle  velocity 
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goes  up  and  down.  It  was  concluded  that  on  average  the  particle  velocity  remains 
constant  with  varying  Al  particle  diameter.  There  is  a  very  small  drop  on  particle 
velocity  from  10%  to  15%  to  20%. 


PBX9501/A1  Particle  Velocity 


-♦-10% 

-•■15% 

^  20% 


Figure  54.  PBX9501  All  Cases  Gas  Particle  Velocity 


Figure  55  plots  the  gas  density  versus  Run  ^  for  the  three  aluminum  mass  frac¬ 
tions.  It  can  be  seen  from  Figure  55  the  gas  density  steadily  increases  as  the  Al 
particle  diameter  increases  at  10%  Al  mass  fraction.  For  15%  Al  mass  fraction  the 
first  density  is  less  than  than  all  the  10%  Al  mass  fraction  cases.  The  same  trend  is 
observed  for  15%  Al  mass  fraction,  where  the  gas  density  steadily  increases  as  the  Al 
particle  diameter  increases.  Again  for  20%  Al  mass  fraction  the  hrst  density  is  lower 
than  all  the  10%  and  15%  Al  mass  fraction  cases.  The  same  trend  is  observed  for 
20%  Al  mass  fraction  cases  where  the  gas  density  steadily  increases  as  the  Al  particle 
diameter  increases. 
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PBX9501/A1  Density 
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Figure  55.  PBX9501  All  Cases  Gas  Density 


4.5  Conclusion 

This  chapter  presented  the  results  from  the  three  main  studies  that  were  per¬ 
formed.  The  purpose  of  these  studies  was  to  evaluate  the  stability  and  limitations 
of  the  MPEXS  code.  The  first  study  was  a  short  study  with  neat  PBX9501  to  com¬ 
pare  run-to-detonation  properties  against  a  single  phase  explosive  code  developed  by 
AFRL.  The  second  study  was  a  convergence  test  of  detonation  properties  by  varying 
the  number  of  computational  cells  for  metalized  PBX9501.  The  third  study  was  to 
study  the  effects  of  aluminum  on  PBX9501  by  performing  several  simulations  at  a 
range  of  aluminum  mass  fractions  and  aluminum  particle  diameters. 
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V.  Conclusions 


To  find  solutions  to  the  bigger  issues  such  as  the  risks  of  injuring  or  killing  inno¬ 
cent  civilians  in  current  and  future  asymmetrical  warfare,  smaller  more  elementary 
problems  need  to  be  understood.  The  main  focus  of  this  research  was  to  study  and 
understand  the  effects  of  aluminum  in  high  explosive  such  as  PBX9501  using  a  com¬ 
putational  simulation  code.  The  MPEXS  hydrocode  was  utilized  in  this  research  to 
perform  simulations  of  a  wide  range  of  metalized  explosive  cases  where  the  size  of  the 
Al  particle  and  the  Al  mass  fraction  was  varied.  The  results  were  separated  into  two 
sections,  the  performance  evaluation  of  the  MPEXS  code  and  the  study  of  the  effects 
of  aluminum  on  PBX9501. 

5.1  Performance  Evaluation  of  the  MPEXS  Code 

The  MPEXS  code  was  evaluated  by  performing  several  studies  to  determine  the 
stability  and  limitations  of  the  code.  The  first  study  was  a  short  study  with  neat 
PBX9501  to  compare  Run2Det  properties  against  the  SPEXS  code  for  different  piston 
speeds.  It  was  important  to  establish  a  baseline  by  studying  the  behavior  of  the  code 
with  neat  PBX9501  with  no  aluminum.  The  second  study  was  a  convergence  test 
of  detonation  properties  by  varying  the  number  of  computational  cells  for  metalized 
PBX9501.  The  goal  of  this  study  was  to  test  the  stability  of  the  code  by  looking  for 
a  convergence  in  the  data  as  the  number  of  computational  cells  increased. 

The  first  study  showed  that  the  average  rate  of  change  for  Run2Det  distance  and 
time  between  each  set  of  mass  fraction  points  was  larger  for  the  SPEXS  code.  The 
difference  in  results  in  both  codes  does  not  imply  the  MPEXS  code  is  bad,  it  only 
means  that  different  equations  of  state  and  burn  rates  were  utilized.  The  takeaway 
is  the  MPEXS  code  followed  the  same  trend  as  the  SPEXS  code,  where  Run2Det 
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distance  and  time  decreases  exponentially  as  the  piston  impact  speed  is  increased. 

The  second  study  showed  that  the  stability  of  the  code  was  proved  by  the  con¬ 
verging  results  as  the  number  of  computational  cells  increased.  The  number  of  com¬ 
putational  cells  is  related  to  the  resolution  or  the  grid  size  of  the  simulation.  It  was 
concluded  that  as  the  resolution  increased,  the  results  became  closer  and  closer  to 
the  same  number,  showing  a  converging  pattern.  The  limits  of  the  MPEXS  were  also 
tested  with  the  small  computational  cell  sizes  of  Nx  =  200  and  Nx  =  400.  When 
these  small  values  were  utilized  the  code  failed  to  run.  The  code  ran  for  Nx  =  800 
and  Nx  =  1000  but  these  small  values  provide  unstable  results  that  varied.  It  was  not 
until  Nx  values  larger  than  1300  when  the  Run2Det  result  began  to  converge.  It  can 
be  concluded  that  the  MPEXS  code,  although  its  limitations,  is  stable  and  provides 
expected  results. 

Below  is  a  short  bulleted  list  to  summarize  this  section  on  performance  evaluation 
of  the  MPEXS  code. 

•  Study  1:  Neat  Explosive 

—  Different  from  SPEXS  due  to  different  EOS 
—  MPEXS  code  follows  same  trend  as  SPEXS 

•  Study  2:  Convergence  Test 

—  Run2Det  distance  and  time  results  converged  as  Nx  increased 
—  Detonation  velocity  results  converged  as  Nx  increased 
—  MPEXS  code  is  stable  and  provides  expected  results 

5.2  The  Effects  of  Aluminum  on  PBX9501  Explosive 

The  main  focus  of  this  research  was  to  use  the  MPEXS  hydrocode  to  run  a  wide 
range  of  simulations  of  metalized  PBX9501  and  determine  what  effects  aluminum 
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particulates  have  on  PBX9501.  Three  aluminum  mass  fractions  were  utilized;  10,  15, 
and  20  percent.  At  each  of  those  aluminum  mass  fractions,  six  aluminum  diameters 
were  tested;  9,  13,  19,  23,  29,  and  34  micron  diameter. 

The  18  simulation  runs  were  performed  for  two  computational  cell,  800  and  2000. 
The  purpose  of  running  simulations  with  both  sizes  of  computational  cells,  800  and 
2000  cells,  was  to  study  the  effects  of  grid  size  on  metalized  explosive  simulations.  It 
can  be  concluded  that  the  simulations  using  2000  computational  cells  provide  more 
consistent  results  for  both  varying  parameters,  Al  mass  fraction  and  Al  particle  di¬ 
ameter.  In  Section  4.3,  it  was  found  that  at  lower  resolutions  the  detonation  property 
results  are  unstable.  The  same  stability  issues  at  the  low  resolution  using  800  cells 
was  seen  especially  for  detonation  velocity  at  15  and  20  percent  Al  mass  fraction. 

As  mentioned  before,  run-to-detonation  properties  are  very  important  in  the  eval¬ 
uation  of  of  explosive  materials.  Run2Det  distance  and  Run2Det  time  were  plotted 
for  all  the  18  simulation  runs.  Increasing  the  Al  mass  fraction  from  10%  to  15%  de¬ 
creased  the  Run2Det  distance  by  an  average  of  0.012  mm,  where  as  increasing  the  Al 
particle  diameter  decreased  the  Run2Det  distance  by  an  average  of  0.001  mm.  The 
same  observation  was  made  for  Run2Det  time.  It  was  concluded  that  varying  the 
aluminum  mass  fraction  had  a  larger  effect  on  Run2Det  properties  than  varying  the 
aluminum  particle  diameter  for  both  Run2Det  distance  and  time.  In  general,  increas¬ 
ing  the  Al  mass  fraction  will  decrease  the  the  time  and  distance  for  the  detonation  to 
reach  steady-state.  If  shorter  Run2Det  distance  and  time  is  a  desired,  then  increase 
the  Al  mass  fraction 

The  other  detonation  properties  such  as  detonation  velocity,  pressure,  and  density 
had  the  following  results  as  Al  mass  fraction  and  particle  diameter  were  varied.  Deto¬ 
nation  velocity  steadily  decreased  as  the  Al  particle  diameter  was  increased.  Also,  as 
the  Al  mass  fraction  increased  the  detonation  velocity  decreased.  The  results  showed 
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that  the  gas  pressure  steadily  increased  as  the  Al  particle  diameter  increased  at  10% 
Al  mass  fraction.  The  same  pressure  increase  with  particle  diameter  was  seen  for 
15%  and  20%  Al  mass  fraction.  For  each  increase  in  mass  fraction  the  first  pressure 
was  less  than  than  all  the  previous  Al  mass  fraction  cases,  clearly  shown  in  Figure 
53.  The  gas  density  showed  the  same  trends  as  gas  pressure  as  discussed  in  this  para¬ 
graph.  The  gas  density  steadily  increased  as  the  Al  particle  diameter  increased  but 
decreased  when  Al  mass  fraction  was  increased.  It  can  be  concluded  that  on  average 
the  particle  velocity  remained  constant  with  varying  Al  particle  diameter.  There  is  a 
very  small  decrease  in  particle  velocity  from  10%  to  15%  to  20%,  but  not  signihcant 
enough  to  vary  the  Al  mass  fraction  in  the  explosive. 

The  addition  of  metal  particulates  in  explosive  mixtures  increases  density  to  the 
shock  wave,  causing  a  higher  pressure  in  the  shock.  The  high  pressure  in  the  shocks 
are  devastating  and  will  incapacitate  adversaries  in  physical  proximity  of  the  detona¬ 
tion  of  the  munition.  The  addition  of  Al  particulates  in  explosive  mixtures  and  low 
fragmentation  weapons  with  concentrated  pressures  that  dissipate  quickly  will  pro¬ 
vide  the  capabilities  needed  to  develop  weapons  needed  to  hght  current  and  future 
asymmetric  warfare.  In  conclusion,  this  research  significantly  contributes  to  both  cur¬ 
rent  explosive  simulation  analysis  and  development  of  future  explosive  formulations 
for  better  munitions  that  will  save  many  lives. 

Below  is  a  short  bulleted  list  to  summarize  this  section  on  the  effects  of  aluminum 
on  PBX9501  explosive. 

•  Aluminum  Particle  Size  -  Increases 

—  Pressure  -  Increases 

—  Density  -  Increases 

—  Particle  Velocity  -  Remained  constant 
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—  Run2Det  Distance  and  Time  -  Small  decrease 
—  Detonation  Velocity  -  Decreases 

•  Aluminum  Mass  Fraction  -  Increases 

—  Pressure  -  Decreases 
—  Density  -  Decreases 
—  Particle  Velocity  -  Remained  constant 
—  Run2Det  Distance  and  Time  -  Decreases 
—  Detonation  Velocity  -  Decreases 

5.3  Recommendations 

Based  on  the  results  for  metalized  explosive  PBX9501,  the  best  explosive  mixture 
is  PBX9501  at  10%  aluminum  mass  fraction  and  34  pm  aluminum  particle  diameter. 
PBX9501_10_34  has  the  highest  gas  pressure  and  density  for  both  particle  diameter 
and  mass  fraction.  The  gas  pressure  for  this  case  is  43.64  GPa  and  the  gas  density  is 
2419  kg/m^.  The  detonation  velocity  was  the  lowest  for  10%  mass  fraction  but  it  was 
higher  than  the  other  mass  fractions.  The  detonation  velocity  for  this  case  is  7.997 
km/ s.  The  Run2Det  time  is  0.258  ps  and  the  Run2Det  distance  is  1.017  pm.  All  these 
properties  are  for  a  piston  speed  of  1.5  km/ s.  If  a  shorter  Run2Det  time  and  distance 
is  desired,  the  aluminum  mass  fraction  can  be  increased  to  15%  with  an  aluminum 
particle  diameter  of  34  pm  but  there  will  be  some  losses  in  other  parameters.  For  15% 
Al  mass  fraction,  the  Run2Det  time  is  0.232  ps  and  the  Run2Det  distance  is  0.926 
pm.  The  pressure  and  density  will  not  affected  as  much  by  the  increase  in  aluminum. 
The  gas  pressure  for  this  case  is  42.31  GPa  and  the  gas  density  is  2361  kg/m^.  The 
detonation  velocity  for  this  case  is  7.848  km/s,  a  loss  of  149  m/s.  The  aluminum 
mass  fraction  of  20%  can  be  utilized  but  the  benehts  gained  are  not  worth  it. 
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Below  is  a  short  bulleted  list  to  summarize  this  paragraph  on  explosive  recom¬ 
mendations. 

•  Best  mixture  -  PBX9501  10%  34  /xm 

—  Highest  Pressure:  P  =  43.64  GPa 

—  Highest  Density:  p  =  2419  kg/m^ 

—  Lowest  Detonation  Velocity  for  10%  but  higher  than  15%  and  20%:  D  = 
7.997  km/s 

—  Run2Det  Time:  t  =  0.258  ps 

—  Run2Det  Distance:  x  =  1.017  mm 

•  If  shorter  Run2Det  time  and  distance  is  desired 

—  Increase  Al  mass  fraction  to  15% 

—  Results  in  losses  in  Pressure,  Density,  and  Detonation  Velocity 

•  No  -  Al  mass  fraction  of  20% 

—  Results  in  losses  in  Pressure,  Density  -  deviates  from  munitions  purpose 

The  MPEXS  hydrocode  is  still  in  a  testing  phase,  therefore  there  are  many  op¬ 
portunities  for  improvement  and  ease  of  use.  Throughout  this  research  there  were 
several  updates  to  the  code  that  came  from  testing.  The  most  current  version  that 
was  utilized  to  generate  the  results  presented  in  this  document  still  has  limitations 
and  issues  to  be  worked  on.  In  Section  4.3  it  was  documented  that  the  code  would 
stop  short  of  the  hnal  run  time  for  all  the  cases.  The  short  running  times  might  be 
due  to  larger  numerical  errors  in  lower  resolution  simulations  versus  errors  in  higher 
resolution  simulations.  The  short  running  time  issue  should  be  corrected  in  future 
versions  of  the  code.  Another  reported  issue  discussed  in  Section  3.6.2  was  that  the 
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code  fails  and  quits  the  simulation  run  for  piston  impact  speeds  of  2.0  km/s  at  high 
Al  mass  fraction  cases.  This  high  metal  mass  fraction,  high  impact  speed  problem  has 
been  reported  to  AFRL,  so  it  can  be  resolved  in  future  versions  of  the  code.  It  would 
be  useful  if  the  MPEXS  code  provided  feedback  to  the  user  when  bad  inputs  are 
entered  or  when  input  parameters  are  outside  the  limits  of  the  code.  Error  messages 
before  the  simulation  begins  should  be  added  to  future  versions  of  the  code.  The 
MPEXS  code  is  stable  and  provides  good  results  for  metalized  PBX9501  explosive. 
The  next  step  is  to  make  the  code  more  user  friendly,  add  more  EOS  to  handle  more 
explosives,  and  improve  feedback  for  inputs  that  will  make  the  code  crash. 

Below  is  a  short  bulleted  list  to  summarize  this  paragraph  on  future  MPEXS  code 
versions. 

•  Fix  all  reported  issues 

—  Simulations  stopping  short  of  hnal  time 

—  Issues  with  high  Al  mass  fraction  with  high  impact  speeds 

•  Make  code  more  user  friendly 

—  Better  input  interface 

—  Error/ Warning  messages  for  inputs  outside  limits 

•  Add  capability  for  more  explosives  and  metals 
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